{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# 1D smooth optimal transport\n",
    "\n",
    "\n",
    "This example illustrates the computation of EMD, Sinkhorn and smooth OT plans\n",
    "and their visualization.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Author: Remi Flamary <remi.flamary@unice.fr>\n",
    "#\n",
    "# License: MIT License\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pylab as pl\n",
    "import ot\n",
    "import ot.plot\n",
    "from ot.datasets import make_1D_gauss as gauss"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate data\n",
    "-------------\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#%% parameters\n",
    "\n",
    "n = 100  # nb bins\n",
    "\n",
    "# bin positions\n",
    "x = np.arange(n, dtype=np.float64)\n",
    "\n",
    "# Gaussian distributions\n",
    "a = gauss(n, m=20, s=5)  # m= mean, s= std\n",
    "b = gauss(n, m=60, s=10)\n",
    "\n",
    "# loss matrix\n",
    "M = ot.dist(x.reshape((n, 1)), x.reshape((n, 1)))\n",
    "M /= M.max()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot distributions and loss matrix\n",
    "----------------------------------\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAADFCAYAAABzYARGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl4U2X2wPHvoWW1yFZEsSCgiGylQEUERHBQdFBRB0cEFVdcBmVUEHQQEZfB0UH9CeIwMo47KCqDK4iAICrSYgsUZEepiiJiEdm6vL8/ToKltCVt09ws5/M8edIkNzcnN2nOfXdxzmGMMcZUVBWvAzDGGBMdLKEYY4wJCksoxhhjgsISijHGmKCwhGKMMSYoLKEYY4wJCksoxhhjgsISijHGmKCwhGKMMSYo4r0OoKjExETXrFkzr8MwxhgDpKen/+ScaxjItmGXUJo1a0ZaWprXYRhjjAFE5OtAt7UqL2OMMUFhCcUYY0xQBJRQRORcEVkrIhtEZHQxj1cXkRm+x5eKSDPf/VVF5HkRWSkia0Tk7uCGb4wxJlwcsQ1FROKAycDZQDawTERmO+dWF9rsOmCnc+4kERkIPAJcBlwKVHfOtReRWsBqEXnVObcl2G/EGFN+ubm5ZGdns2/fPq9DMR6pUaMGSUlJVK1atdz7CKRRvguwwTm3CUBEpgP9gcIJpT8wzvf3TGCSiAjggKNEJB6oCRwAdpU7WhOQOXNg9WoYPhyqWKWmCUB2dja1a9emWbNm6L+uiSXOOXbs2EF2djbNmzcv934C+bk5Htha6Ha2775it3HO5QE5QAM0ufwGfA98AzzmnPu56AuIyFARSRORtO3bt5f5TZjfPfMM/PGPcMcdMGgQ7N/vdUQmEuzbt48GDRpYMolRIkKDBg0qXEKt7PPXLkA+0BhoDtwpIi2KbuScm+qcS3XOpTZsGFB3Z1OEczBmDNx8M5x3Hjz4IMyYoX/n5HgdnYkElkxiWzA+/0CqvL4FmhS6neS7r7htsn3VW3WAHcAg4APnXC7wo4gsAVKBTRUN3Bzqtttg0iS49lr4178gPh6aNtXbPXvCkiWQkOB1lMaYaBZICWUZ0FJEmotINWAgMLvINrOBIb6/BwDznS5W/w1wFoCIHAV0Bb4KRuDmd19/DZMnw9Ch8OyzmkwArrwSXn8dVqyAV17xNkZjjuShhx6ibdu2JCcnk5KSwtKlS70O6TAJvrOy7777jgEDBpS43S+//MLTTz9d6r66desGwMKFCzn//PPLFMesWbNYvfr3ZuyxY8cyb968Mu2jMhwxofjaRIYBc4A1wGvOuSwRGS8iF/o2mwY0EJENwB2Av2vxZCBBRLLQxPScc25FsN9ErHvuOb2++24oWmrt3x/atoVp00IflzGB+uyzz3jnnXdYvnw5K1asYN68eTRp0uTITzyCvLy8IER3uMaNGzNz5swSHy8tofhj+vTTT8v9+kUTyvjx4+nTp0+59xcsAU294px7D3ivyH1jC/29D+0iXPR5u4u73wRPfr4mlD59oLgp0ETguuu0kX7VKmjXLuQhmgjz179CRkZw95mSAk88UfLj33//PYmJiVSvXh2AxMTEg4999NFHjBgxgry8PE499VSmTJlC9erVD07TlJiYSFpaGiNGjGDhwoWMGzeOjRs3smnTJpo2bcpLL73EqFGj+OCDD6hSpQo33HADt956K+np6dxxxx3s3r2bxMRE/vvf/3LccccdEtfmzZsZNGgQu3fvpn///gfv37JlC+effz6rVq0iKyuLa665hgMHDlBQUMAbb7zBvffey8aNG0lJSeHss8+mX79+3HvvvdSrV4+vvvqKdevWkZCQwO7duwHYtWsX/fr1Y8OGDfTu3Zunn36aKlWqHLLNzJkzeeeddxg6dCizZ8/m448/5sEHH+SNN97ggQce4Pzzz2fAgAGlHq8hQ4bw9ttvk5uby+uvv84pp5wSrI8YsJHyEe+jj+CbbzRplOTKK6FqVSulmPB1zjnnsHXrVk4++WRuueUWPv74Y0B7n1199dXMmDGDlStXkpeXx5QpU464v9WrVzNv3jxeffVVpk6dypYtW8jIyGDFihUMHjyY3Nxcbr31VmbOnEl6ejrXXnstf/vb3w7bz/Dhw7n55ptZuXLlYcnG75lnnmH48OFkZGSQlpZGUlISEyZM4MQTTyQjI4NHH30UgOXLl/Pkk0+ybt26w/bxxRdf8NRTT7F69Wo2btzIm2++WeJ769atGxdeeCGPPvooGRkZnHjiiQcfO9LxSkxMZPny5dx888089thjRzyOZRV2k0Oaspk2DerXh4suKnmbxESt+nrxRZgwAXwngcYUq7SSRGVJSEggPT2dxYsXs2DBAi677DImTJhAx44dad68OSeffDIAQ4YMYfLkyfz1r38tdX8XXnghNWvWBGDevHncdNNNxPsaF+vXr8+qVatYtWoVZ599NgD5+fnFJowlS5bwxhtvAHDllVcyatSow7Y5/fTTeeihh8jOzuaSSy6hZcuWxcbUpUuXEsd4dOnShRYttAPs5ZdfzieffFJqG01J1q5dW+rxuuSSSwDo3LlzqUmrvKyEEsF27IBZs7QEcqQkcf31uv3sot0pjAkTcXFx9OrVi/vvv59JkyYd/CEvSXx8PAUFBQCHjZ846qijSn2uc462bduSkZFBRkYGK1euZO7cucVue6TutIMGDWL27NnUrFmTP/7xj8yfP7/Y7UqLqehr+G8Xvj8Ysxj4qxTj4uIqpX3JEkoEe+klOHCg9Oouvz59oEkTq/Yy4Wnt2rWsX7/+4O2MjAxOOOEEWrVqxZYtW9iwYQMAL774ImeeeSagS12kp6cDlJp8zj77bP71r38d/AH9+eefadWqFdu3b+ezzz4DdOqZrKysw57bvXt3pk+fDsDLL79c7P43bdpEixYtuO222+jfvz8rVqygdu3a/PrrrwG//y+++ILNmzdTUFDAjBkz6NGjBwCNGjVizZo1FBQU8NZbbx3cvqT9l3a8QsESSoRyTrsIn3oqtG9/5O3j4uCaa2DuXG1zMSac7N69myFDhtCmTRuSk5NZvXo148aNo0aNGjz33HNceumltG/fnipVqnDTTTcBcN999zF8+HBSU1OJi4srcd/XX389TZs2JTk5mQ4dOvDKK69QrVo1Zs6cyahRo+jQoQMpKSnF9rp68sknmTx5Mu3bt+fbb4sOv1OvvfYa7dq1IyUlhVWrVnHVVVfRoEEDunfvTrt27Rg5cuQR3/+pp57KsGHDaN26Nc2bN+fiiy8GYMKECZx//vl069btkCq5gQMH8uijj9KxY0c2btx48P7SjlcoiA4XCR+pqanOFtg6spUrITkZpkyBQL8vW7ZA8+bajlJMVbCJYWvWrKF169Zeh2E8Vtz3QETSnXOpgTzfSigRyl9N269f4M9p1kzHpCxYUCkhGWNinCWUCLVgAZx4oraLlEXv3vDJJ5CbWzlxGWNilyWUCFRQAIsWQa9eZX9ur17w229gtYrGmGCzhBKBMjNh587yJRR/hw+r9jLGBJsllAi0cKFelyehJCZqrzD/PowxJlgsoUSgBQvgpJMgKal8z+/VS6ezP3AgqGEZY2KcJZQIk5+v7Se9e5d/H717w549sGxZ8OIypiJ27NhBSkoKKSkpHHvssRx//PEHbx+opDOf5cuX88EHHwS0bY8ePcjwzZjZt2/fUgctTpw4sdRR7ddccw1r164lLy+PunXrVijmt9566+BcYeHA5vKKMBkZugJjeaq7/Hr21OsFC6B796CEZUyFNGjQ4OAP9rhx40hISGDEiBEBPz8/P7/UwY3FWb58OatWreLcc88t0/PmzJlT6uMTJ07k2muvpUaNGsXG+ZxvvYnyTH1SNGb/AMhwYQklwlSk/cSvQQMdFLlwoS4bbMwhvJi/vhQXXHAB3333Hfv27eP222/n+uuvJy8vj8TERK6++mrmz5/Pv/71L7Zv387IkSNJSEigW7dubN26lVmzZrF7926GDRvG6tWryc3NPbh2yPjx49m7dy8LFy5kzJgxh0zGuGfPHoYMGcKqVato06bNISWOpKQkVq1aRVxcHH/+85/57rvvyM/PZ9y4cWzdupUff/yRM844g0aNGvHBBx8cFufIkSOZNGkS7XxrSdx222189NFHNG7cmOnTp9OgQQN69OjBpEmTSElJYdu2bfTo0YOVK1ceFvMvv/zCqlWreOKJJ9i8eTPXXnstO3bsoFGjRjz33HMkJSVxxRVX0KBBA5YtW8a2bdv45z//WWmJyKq8IsyCBXDyydC4ccX207u3tqPs3x+cuIypLM8//zzp6eksW7aMiRMnsnPnTgBycnLo2bMnK1asoEOHDtxyyy3MnTuXtLQ0tm3bdvD548eP59xzz+WLL75g/vz53HnnnYgIY8eOZfDgwWRkZBw2s++kSZOoV68ea9asYcyYMXz55ZeHxfXee+/RrFkzMjMzD85cfPvtt3PMMcewePHigysoFo7z9NNPP2QfOTk5dO/enaysLE4//XQeeOCBEo9DzZo1S435lltu4frrr2fFihVceumlh8zI/OOPP7JkyRJmzZrF3XffHeCRLzsroUSQvDxYvBgGDqz4vnr1giefhC++gDPOqPj+TBTxYv76Ujz++OPM9k2TnZ2dfXDhqmrVqh080169ejWtWrXihBNOAHQK+BdeeAGAuXPn8v777zNhwgRAZ+395ggT2i1atIi77roLgI4dO9K2bdvDtklOTmb06NGMHj2aCy64gO4l1B8XjrOo+Ph4Lr1U1yC84oorGDRoUKlxlWbp0qW88847AFx11VXce++9Bx+76KKLEBGSk5NLnJMsGKyEEkG+/BJ27apYg7xfz566mqONRzHhbN68eSxatIjPP/+czMxMkpOTD1Y/1axZ84hTy4NOVT9r1qyDU9V/8803B9cLqYjWrVuTlpZG27ZtGT16NA8//HCx2wUaJ/w+XX1pU/OXR/VC61tU5vyNllAiyKJFeh2M2ajr19d2FP8+jQlHOTk51K9fn5o1a5KVlcWyEromtmnThrVr17J161acc8yYMePgY3379uWpp546eNtffVXaFPM9e/bklVdeASAzM7PYqe2//fZbEhISuPLKK7nzzjtZvnz5EfdbVF5e3sGFrl555ZWD09YXnpq/8Nr1pe27a9euvPbaawC89NJL9PT3vgkhSygRJC0NmjaFElYiLbOuXXWfYTbhtDEH9evXjz179tCmTRvGjBnDaaedVux2tWrVYtKkSfTp04fU1FTq1q1LnTp1AJ3m/rfffqN9+/a0bduWcePGAXDWWWeRmZlJx44dD/nRBhg2bBg7duygdevWPPDAA3Ts2PGw18zMzOTUU08lJSWFhx9+mHvuuQeAoUOH0qdPH/r06XPE91enTh0WL15M27Zt+eSTTxjj6yUzcuRInnzySTp16nSwzehIMU+ePJmpU6eSnJzMjBkzePzxx4/4+sFm09dHkJYttVRxhIXsAvbss3DDDbB+vQ6UNLErGqav3717NwkJCTjnuPHGG2nfvj233nqr12FFFJu+Pkbs3AkbNkBqQB9rYPz7svxtosGUKVNISUmhTZs27N27lxtuuMHrkGKO9fKKEL7q2aAmlLZtdS36tLTg9BwzxksjR44MaHVEU3mshBIh/KWIzp2Dt8+qVXW8mZVQDFRu7x8T/oLx+VtCiRBpadCihfbOCqbUVEhP1zVWTOyqUaMGO3bssKQSo5xz7Nixo9jpYsrCqrwiRFoadOkS/P2mpsLkybBuHZxySvD3byJDUlIS2dnZbN++3etQjEdq1KhBUnmnMPexhBIBfvoJtmyBm28O/r79VWhpaZZQYlnVqlVp3ry512GYCGdVXhHAN74pqA3yfq1bQ82a1o5ijKm4gBKKiJwrImtFZIOIjC7m8eoiMsP3+FIRaVbosWQR+UxEskRkpYhUrJIuBvkTSqdOwd93fDx07Pj7axhjTHkdMaGISBwwGTgPaANcLiJtimx2HbDTOXcS8DjwiO+58cBLwE3OubZALyA3aNHHiLQ0HdRYxrV4Apaaqt2S8/MrZ//GmNgQSAmlC7DBObfJOXcAmA70L7JNf+B5398zgT+IznJ2DrDCOZcJ4Jzb4Zyzn60ySkurnOouv9RUXcHxq68q7zWMMdEvkIRyPLC10O1s333FbuOcywNygAbAyYATkTkislxE7iruBURkqIikiUia9TI51A8/wNatlZ9QwNpRjDEVU9mN8vFAD2Cw7/piEflD0Y2cc1Odc6nOudSGDRtWckiRpTIb5P1OPhkSEiyhGGMqJpCE8i3QpNDtJN99xW7jazepA+xASzOLnHM/Oef2AO8BldC0HL3S0nTdkmImOw2auDht8LeEYoypiEASyjKgpYg0F5FqwEBgdpFtZgNDfH8PAOY7HXI7B2gvIrV8ieZMYHVwQo8NaWnQqhXUrl25r5OaqsuI5+VV7usYY6LXEROKr01kGJoc1gCvOeeyRGS8iFzo22wa0EBENgB3AKN9z90JTESTUgaw3Dn3bvDfRvTKyKic7sJFdewI+/bB2rWV/1rGmOgU0Eh559x7aHVV4fvGFvp7H3BpCc99Ce06bMro55+1Qb5Dh8p/Lf9rZGbqLMTGGFNWNlI+jK1YodehSCinnALVqmlCMcaY8rCEEsYyMvQ6JaXyX6tqVS2Z+F/TGGPKyhJKGMvMhEaN9BIKHTpYCcUYU36WUMJYZmZoqrv8OnTQgZTbtoXuNY0x0cMSSpjKzYWsrNBUd/n5X8tKKcaY8rCEEqa++goOHAh9CQUsoRhjyscSSpjy/6iHMqHUqwdNmlhCMcaUjyWUMJWRAdWr6yj5UEpJsZ5expjysYQSpjIzoV07XQArlDp00NHy+/aF9nWNMZHPEkoYci70Pbz8OnTQhbayskL/2saYyGYJJQxt2wbbt4e2h5ef9fQyxpSXJZQw5G/D8KKE0qKFro1i7SjGmLKyhBKG/KWD5OTQv3aVKtC+vZVQjDFlZwklDGVmQrNmULeuN6+fkqIxOOfN6xtjIpMllDCUkeFNdZdfhw6QkwNff+1dDMaYyGMJJczs3Qvr1nmfUMCqvYwxZWMJJcysWgUFBd708PJr317XsbeEYowpC0soYcbLHl5+Rx0FLVtaTy9jTNlYQgkzmZlQu7Y2ynvJ1kYxxpSVJZQwk5mp3YWrePzJdOgAmzbBrl3exmGMiRyWUMJIQYEmFC/bT/z8MfjXtTfGmCOxhBJGtmyBX3/1tv3Ez3p6GWPKyhJKGPFiDZSSHH881K9vCcUYEzhLKGEkM1PbTtq18zoS7TZsDfPGmLKwhBJGMjK0u26tWl5Hojp0gJUrdTp7Y4w5EksoYSRcGuT9UlJ05P769V5HYoyJBJZQwsQvv2ijfDi0n/hZw7wxpiwsoYQJf/fccEoorVvrEsQ2Yt4YE4iAEoqInCsia0Vkg4iMLubx6iIyw/f4UhFpVuTxpiKyW0RGBCfs6BNOPbz8qlfXpGIlFGNMII6YUEQkDpgMnAe0AS4XkTZFNrsO2OmcOwl4HHikyOMTgfcrHm70ysyExERo3NjrSA7lXxvFGGOOJD6AbboAG5xzmwBEZDrQH1hdaJv+wDjf3zOBSSIizjknIhcBm4HfghZ1FPKvgSLidSSH6tABXnxR17hv2NDraExA9u2DtWth61a95OTomUqTJnDCCdC8efh90UxUCCShHA9sLXQ7GzitpG2cc3kikgM0EJF9wCjgbKDE6i4RGQoMBWjatGnAwUeLvDydtv4vf/E6ksMVbpjv08fbWEwp9u6FDz6A11+Ht9+G3btL3rZFC/jzn/WSkmLJxQRNZTfKjwMed86V8u0G59xU51yqcy61YQyeBq9bB/v3h1f7iZ/19Apz+/fD449DUhJccgnMnQuDBsGMGfD555Cdrcll3Tr46COYMkUHOz36KHTqBF27wqJFXr8LEyUCKaF8CzQpdDvJd19x22SLSDxQB9iBlmQGiMg/gLpAgYjsc85NqnDkUeTLL/U6nMag+DVsqLUl/hhNmHAOpk+He+7R/ubnnAN33glnnaVd84pq2VIvZ50FN90EP/0Er70GDz8MZ54JF1wAjzyivTCMKadASijLgJYi0lxEqgEDgdlFtpkNDPH9PQCY79QZzrlmzrlmwBPAw5ZMDpeeDjVqQJuiXR3CROfOGqMJEz//DBddpCWROnVgzhy9nHNO8cmkOImJcMstWnL5+9/h44/1jObppzVZGVMOR0wozrk8YBgwB1gDvOacyxKR8SJyoW+zaWibyQbgDuCwrsWmZOnpWrUU6G9BqHXurG28v/7qdSSGzz6Djh3h/fdh4kRYvlwTSXnVqgWjR+t0CH36aEPepZfqSFtjyiignzDn3HvAe0XuG1vo733ApUfYx7hyxBf1Cgq0OunKK72OpGSdO+tJa2Ym9OjhdTQxbMoUuO027a21ZAmcemrw9n3MMdqYP3Ei3H23Jqp337UqMFMmNlLeYxs26Jl/p05eR1Iyf2xW7eUR52DcOK2iOu88PQMJZjLxq1IFRoyAxYu119gZZ8AXXwT/dUzUsoTiMf+PdOfO3sZRmsaN4dhjLaF4Ij8fhg2D+++Ha6+FN9/UdpPK1LWrloDq1NFG/A8/rNzXM1HDEorH0tN1ipO2bb2OpHTWMO+B/Hy46iptKL/rLnj22dA1tLVooUnlpJOgXz+YNSs0r2simiUUjy1fDsnJULWq15GUrnNn+Oor+M3mOwgN5+Dmm+GVV7Rr7yOPhH4A4rHHwsKF+uFfdhnMmxfa1zcRxxKKh5zThBLO1V1+nTtrBwIb4BgCzmmJ5N//hr/9TRvJvVK3Lrz3HrRqpV2VP/vMu1hM2LOE4qGNG3WapUhJKGDVXiHx97/DY49pF94HHvA6GqhXT0fgH3cc/PGPdlZhSmQJxUOR0CDv17gxNGpkCaXSvfCClkquuAL+7//CZ56tY4/VKq+EBE0q3xadLMMYSyieSk+HatXCv0Ee9HetUydLKJXqk0/g+uu1Z9V//qPdeMPJCSdo9deuXXDhhdagZg4TZt/Y2LJ8ObRvr0klEnTuDKtXw549XkcShTZtgosv1qnlZ84M314a7dvrHGIZGdoDraDA64hMGLGE4pFIapD38zfM+5crNkGSk6OTM+bnwzvvaJtFOOvXD/75Tx0TM2aM19GYMGIJxSObN8POnZGXUMCqvYKqoEDn3Vm3Dt54Q2cEjgTDh8PQodqB4PXXvY7GhAlLKB6JpAZ5v6Qknc7eEkoQTZjw+xxavXt7HU3gROCpp+D00+Gaa2DNGq8jMmHAEopH0tO1mrxdO68jCZyIjZgPqrlztcpo0CCdXiXSVKumpZOjjtL2n127vI7IeMwSikf8s5BXr+51JGVz2mm6XLFNZV9BW7bA5ZdrF7+pU8One3BZHX+8rg65YYOWVGwtlZhmCcUDubk6iWu3bl5HUnbdumm1/9KlXkcSwQ4c0PXc8/K0Yfuoo7yOqGJ69dKpYd58U5cjNjHLEooHMjJg3z7o3t3rSMqua1c9mf70U68jiWCjRsGyZfDcc5HTCH8kd9yh1V6jRtnZRgyzhOIB/49xJJZQjj5ahyJYQimnWbPgiSd0oaxLLvE6muARgWnTtOfGZZdpF0YTcyyheODTT3XQcePGXkdSPt26aRuQjWkroy1btJ2hc2f4xz+8jib46tXT9pTvvrP2lBhlCSXEnNNlJiKxdOLXrZt26MnK8jqSCJKbCwMHahZ+7bXI640RqC5dNFn+7386F5mJKZZQQmzrVp1XL9ITCli1V5mMGaNtC9Om6eJV0Wz4cJ3r6667dLliEzMsoYSY/0c4Ehvk/Vq00JmHLaEEaO5cPWu/8UYYMMDraCqfiE5uecwx2p5ifcxjhiWUEPv0U+0l2r6915GUn4iWUiyhBGDbNp1apV272OpS26ABvPyyLvoTiYM2TblYQgmxJUt0cGColgavLN266Vi2H37wOpIwVlCgM/L++qvO0FuzptcRhVbPnjB2rK7x8uKLXkdjQsASSgjt3q2L3UVy+4mf/z3YirCleOwx+PBD7SYcCYveVIYxY+DMM+GWW2D9eq+jMZXMEkoILVumM5RHQ0Lp1EmncrJqrxIsXaorL156Kdxwg9fReCcuDl56Sb8sAwfC/v1eR2QqkSWUEPL/+J5+urdxBEONGpCaagmlWL/8oj+exx8f2fN0BUtSks4KsHw53H2319GYSmQJJYQ+/VRrPurW9TqS4OjWDdLS7KTzEM5pb66tW+HVV6Pnw66oCy+EW2/Vjgnvvut1NKaSWEIJkdxcWLwYevTwOpLg6dFDk8nnn3sdSRiZNk0HLj74YHQURYPpH/+AlBQYMkQHY5moE1BCEZFzRWStiGwQkdHFPF5dRGb4Hl8qIs18958tIukistJ3fVZww48cn3+unX3OOcfrSIKnd2/trTZnjteRhImVK/UsvE8fHdRnDlWjhvZ227dP14DJy/M6IhNkR0woIhIHTAbOA9oAl4tImyKbXQfsdM6dBDwOPOK7/yfgAudce2AIELN9B+fM0fbJP/zB60iC5+ij9STcEgrahe/Pf9YqrpdegipW+C9Wq1bwzDOwaBHcf7/X0ZggC+Rb3wXY4Jzb5Jw7AEwH+hfZpj/wvO/vmcAfREScc186577z3Z8F1BSRKJ3EqHRz5+rU73XqeB1JcPXtq22t27d7HYnH/vIXWLtWB/M1auR1NOHtiivg2mvhoYdg3jyvozFBFEhCOR7YWuh2tu++YrdxzuUBOUCDItv8CVjunDusCVdEhopImoikbY/CX6afftLG62iq7vLzv6cPP/Q2Dk/99786eG/sWDgrZmt1y+app6BNGxg8GL7/3utoTJCEpFwuIm3RarAbi3vcOTfVOZfqnEtt2LBhKEIKqXnztPNP375eRxJ8nTrpLBsxW+21YoUO2uvVC+691+toIketWtp5Yfdu7WJt7SlRIZCE8i3QpNDtJN99xW4jIvFAHWCH73YS8BZwlXNuY0UDjkRz5+pSEampXkcSfHFx2gY9d24MLn+RkwN/+pO2m7z6qh4ME7g2bXSczqJFcM89XkdjgiCQhLIMaCkizUWkGjAQmF1km9nZ9NAgAAAQpklEQVRoozvAAGC+c86JSF3gXWC0c25JsIKOJM7p2XufPtH7e9O3r86BuHKl15GEkHNw9dWwebOeaR97rNcRRabBg7WE9+ijuia9iWhHTCi+NpFhwBxgDfCacy5LRMaLyIW+zaYBDURkA3AH4O9aPAw4CRgrIhm+yzFBfxdhLCtLF7CLxuouP387SkxVez32mC7n++ij0TW4yAsTJ+rCXFdfDevWeR2NqQBxYVZPkZqa6tLS0rwOI2j++U8YMQK++QaaNDny9pGqXTs9SY+JTjvz5sG558LFF2vpJNanVgmGb77RBrlGjXTQVu3aXkdkfEQk3TkXUIW9dZavZHPnQuvW0Z1MQEtgixfDnj1eR1LJNm7U8SannKKLSFkyCY6mTTU5r12r68cUFHgdkSkHSyiVaO9ebW+M5uouv7594cABWLjQ60gq0a+/Qv/+mkRmz7az6GA76yyd6+t//4P77vM6GlMOllAq0bvv6iwT/fp5HUnlO+MM/X194w2vI6kkBQV65vzVV3omHe3rwntl2DAd9Pjgg/D6615HY8rIEkolevllbVfo3dvrSCpfzZrapDBzpibRqHPPPXrmPHFidM2fE25E4OmndSrrIUPgiy+8jsiUgSWUSrJzJ7z3no7ZitbuwkUNHgy7dun7jirPPAOPPAI33aSTP5rKVb06vPWWno1dcIF2zTYRwRJKJXnjDW1TGDTI60hC56yz4JhjtGQWNd55R+fp6tdPpwuxRvjQOOYYeP99XffhvPPg55+9jsgEwBJKJXnlFWjZMjpHx5ckPl5LZO++q4sWRrz0dLjsMl3DY/p0fYMmdFq10mrGzZu1M0RU1qVGF0soleDbb7W306BBsXdCO2iQLroV8YOeV6/WsSYNG2opJSHB64hi0xlnwPPPwyefaHft3FyvIzKlsIRSCaZP15k5Yqm6y69LFzjxRC2hRaxNm+Dss7VE8uGHcNxxXkcU2wYOhMmT4e234aqrID/f64hMCSyhVIKXX9aqrpNP9jqS0BPRRDp/vk45E3Gys7UX1759mkxatvQ6IgM639cjj+jZ2o032sDHMGUJJcjWrIEvv9QeT7Fq8GAtoc2Y4XUkZeRPJjt26MRk7dp5HZEp7K67YMwYmDZNx6tYUgk7llCC7JlntKbkssu8jsQ7rVrBqafqsYiY2olNm7S+/vvvtXdRLPWmiCTjx2timTIFrrsugr5gscESShD9+CP8+986oDrWq91HjtSJY996y+tIAvDVV9Czpw6imT8funf3OiJTEhGYMAHGjdOVMgcPtob6MGIJJYiefFKr3keN8joS711yibYhPfxwmC+8lZYGZ56pP0oLF1rJJBKI6Fxfjz6q9aoXX6wrPxrPWUIJkpwcmDQJBgzQKp9YFxcHo0dre1LYrpPy1ltaMqlVS2fxbN/e64hMWYwYofWq77+vn+O3RReSNaFmCSVInn5aa0zuvtvrSMLH4MGQlKSllLDinC6Q9ac/QXKyrr9hZwGR6cYbdZzQ+vVw2mmQkeF1RDHNEkoQ7Nmjs26fey507Oh1NOGjWjVtS1m8WC9h4bffdGXAkSO1OLlggS7qZCLXeefpwEcRXT0zogdBRTZLKEHw7LOwfbtOSGsOdf31kJgIDz3kdSTo6PcuXeDFF7VRd/p0nSbZRL4OHXRm4o4dtWh84402VYsHLKFU0Lffavtgr17a69QcqlYtLQzMmePhWinOwXPPaV/mn37SZTTvuw+q2Nc/qhx3nJY4R42CqVOha1cdGGZCxv6jKsA5uOEGnbtq6lSvowlft9+uy4XffLOW5EJq61adKfjaazWhfPkl9OkT4iBMyMTHa7fid97RgaopKfD3v0NenteRxQRLKBXwn/9oB5MJE2yGjtJUrarz++XkaFIJSTfiggLN8m3bwscfa5/u+fOhceMQvLjxXL9+kJWl66ncc4+WVqzBvtJZQimnb77RM+9evXQWCFO6du3g/vu12qvSp2T5+GMdT3LjjXq9ciXcdptVccWaRo10CdHXX9d/2E6dtErhhx+8jixq2X9YORw4oB2FnNNSiv1OBWbECO3Z+Ze/wJYtlfACWVnaFbhXL20reeUV+OgjW/891g0YAGvXwvDhOrr+pJO0l8iuXV5HFnXsp7CMDhzQZRkWLNAF/Jo39zqiyBEfr1VfBQXQu3cQk0p6ug7Nb9dOG9wffFB/QC6/PPYWpDHFq1dP+/ZnZenSomPGwAknaOcMWw0yaCyhlIE/mfzvf5pMrr7a64giT6tWOiv8L79UMKkcOKB1Z717a7XWggUwdqzu8G9/s+7Apngnn6z/wMuW6Xdn/Hho2hRuukk7bJgKsYQSoD17Dk0m1m5SfqmpMG+eJpVevWDDhgCf6JzOvTVyJDRpogsvbdmi62R8/bU20jRoUImRm6iRmqrLiq5cqf/YL7ygbSynnabTXlg7S7mIC7OZ+1JTU11aWprXYRxi3jwYOlSXtp40SdsATMWlp+vCiPv3wwMPaLv5Ycu2HzgAn36q3elef10/hPh47cVz8826A2vEMhW1c6cOeJ06VavFqlTRSUMvvhj69tVunDFafSoi6c65gGZNtYRSim3bdG6u//5XS8r//rfOQWeCJztbF+N7+209aZzy+D5SJV2TyOLFWpW1e7cmkT/8Qc8mL7oI6tf3OnQTjZzThPL661qlunat3t+smZ68dOumlxhKMEFPKCJyLvAkEAc865ybUOTx6sALQGdgB3CZc26L77G7geuAfOA251ypc896nVDy87Vd99//1h8553Q9n7FjoUYNz8KKPvn52pVz/Xrcqiy2zF7BriUraJ23kmro+hYFLU6iSt+z4ZxztCH16KM9DtrEnI0b9Qdhzhxd3iAnR+9v0EAHTSYn6yzVrVtrkqlfP+oSTVATiojEAeuAs4FsYBlwuXNudaFtbgGSnXM3ichA4GLn3GUi0gZ4FegCNAbmASc750pcZi2UCaWgQBfo27BBJ5xdskRPjHfsgIYNYcgQreqyQYtlsG+fVh/88ot23f3xR718/70WR7Zu1USyefOhCyMdeyy5bZLJlI68sP50ZnzTlZzqjTj1VF3vqnt3/Z9t2lQnnTQm5AoKdCqXzz6DpUshMxNWrYK9e3/fpm5d7frZpIlekpJ0PMwxx+ilfn3dpm5dXeMhAgQ7oZwOjHPO9fXdvhvAOff3QtvM8W3zmYjEA9uAhsDowtsW3q6k16tIQtm+9mcyH36X/HwOXnJz4UAuHNgP+/bDb7t1wtldu/T3LrfQjAzHHavJI7kDdO5UTH1+ZSrtcyj82JH+Lnxd9FJQ8Pt14UvhA5afr9NU5OXpwfNf9u/Xy759etm7V3sq7NkDv/76+2X//uLfQ5UqOtdSUpL+o514oh7sE0+ENm30n63QW1m6VGsdliyB5ct/zz1VqugujjtO/zfr14c6dbRTl/9Stape4uP1fzYuTp8ncvjFr6S/jSmNFOST8MNGjt62jto/bODoH9aTsH0ztXZmU2vHVqrv+aXE5+bWSCC3Rm1ya9Qmr0YCedVqkV+tJvlVa5JftQb5VatTULUG+fHVcHFVKfBdXJU4CuLicVXicBKn11XicFIFqlTBSRWcCPivq1al69Tryv8ey5BQAvnJPB7YWuh2NnBaSds45/JEJAdo4Lv/8yLPPb6YgIcCQwGaNm0aSNzF+in9a/q8cFW5n8823yVcplr3ioj+GletqsWBqlWhevVDL7VqQe3amghq1/79Ureu9vmvV0+nGfafmSUmBpyhRXSmjK5d9fbevZpU1q/Xgs3mzdoJZ/t2reLOydFt9u4N89UhTRSKA072XQ5Xi99oyHaO4Uca8QP12HnwcvS+XdTe9ytHs4sEdlOTvdTkV2rxA9XZT3X2U5P9VOMAVck9eB1PiRU8xdpDTahAQimLUJ6Dl8g5NxWYClpCKe9+TrqwDT99voH4eA5eqlWLoE5ApZ0aB3I67f+78LX/UvgUvfApu/9v/6m8/3YYqVnz92qv0jinncJyc38vYOXnH1oYK1xgK/y84v42puKO8l2alfmZBcBe3+UQ/lqG/HykIP/368L383uthAicUNG3EaBAEsq3QJNCt5N89xW3TbavyqsO2jgfyHODpmpCdRJPO7Gydm/CnMjvBShjopegJaPwa4MJ5FR0GdBSRJqLSDVgIDC7yDazgSG+vwcA8502zswGBopIdRFpDrQEvghO6MYYY8LJEUsovjaRYcAcNCX+xzmXJSLjgTTn3GxgGvCiiGwAfkaTDr7tXgNWA3nAX0rr4WWMMSZy2cBGY4wxJSpLL6/wan01xhgTsSyhGGOMCYqwq/ISke3A1xXcTSLwUxDCiWR2DJQdBzsGfnYcyncMTnDONQxkw7BLKMEgImmB1vlFKzsGyo6DHQM/Ow6VfwysyssYY0xQWEIxxhgTFNGaUKZ6HUAYsGOg7DjYMfCz41DJxyAq21CMMcaEXrSWUIwxxoSYJRRjjDFBEVUJRUTOFZG1IrJBREZ7HU+oiEgTEVkgIqtFJEtEhvvury8iH4rIet91Pa9jrWwiEiciX4rIO77bzUVkqe87McM3wWlUE5G6IjJTRL4SkTUicnqsfRdE5Hbf/8IqEXlVRGrEwndBRP4jIj+KyKpC9xX72Yv6P9/xWCEinSr6+lGTUHxLFU8GzgPaAJf7liCOBXnAnc65NkBX4C++9z4a+Mg51xL4yHc72g0H1hS6/QjwuHPuJGAnEJqVhrz1JPCBc+4UoAN6PGLmuyAixwO3AanOuXbopLYDiY3vwn+Bc4vcV9Jnfx46A3xLdIHDKRV98ahJKOi69Rucc5uccweA6UB/j2MKCefc98655b6/f0V/QI5H3//zvs2eBy7yJsLQEJEkoB/wrO+2AGcBM32bxMIxqAP0RGcAxzl3wDn3CzH2XUBnUq/pW5+pFvA9MfBdcM4tQmd8L6ykz74/8IJTnwN1ReS4irx+NCWU4pYqPmy54WgnIs2AjsBSoJFz7nvfQ9uARh6FFSpPAHehi92BLkP9i3Muz3c7Fr4TzYHtwHO+qr9nReQoYui74Jz7FngM+AZNJDlAOrH3XfAr6bMP+m9mNCWUmCciCcAbwF+dc7sKP+Zb8Cxq+4iLyPnAj865dK9j8Vg80AmY4pzrCPxGkeqtGPgu1EPPvpsDjdE1eItWA8Wkyv7soymhhHS54XAjIlXRZPKyc+5N390/+IuwvusfvYovBLoDF4rIFrS68yy0LaGur9oDYuM7kQ1kO+eW+m7PRBNMLH0X+gCbnXPbnXO5wJvo9yPWvgt+JX32Qf/NjKaEEshSxVHJ11YwDVjjnJtY6KHCSzMPAf4X6thCxTl3t3MuyTnXDP3s5zvnBgML0GWpIcqPAYBzbhuwVURa+e76A7piasx8F9Cqrq4iUsv3v+E/BjH1XSikpM9+NnCVr7dXVyCnUNVYuUTVSHkR+SNaj+5fqvghj0MKCRHpASwGVvJ7+8E9aDvKa0BTdEmAPzvnijbYRR0R6QWMcM6dLyIt0BJLfeBL4Arn3H4v46tsIpKCdkyoBmwCrkFPHmPmuyAi9wOXoT0gvwSuR9sHovq7ICKvAr3Qaep/AO4DZlHMZ+9LtpPQ6sA9wDXOuQotlxtVCcUYY4x3oqnKyxhjjIcsoRhjjAkKSyjGGGOCwhKKMcaYoLCEYowxJigsoRhjjAkKSyjGGGOC4v8B1ZcTbJpyYBQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 460.8x216 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFgCAYAAACFYaNMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XmcHHWd//HXp3tyEBIIEAyBAANyLSSAIcrNAkFBwI0/UeRQoiL3fSyni6DLobIQlgU0giy3CAgIcqywIAgkmhCW+wwEEo4EDHckycz398enerpn0pPpmenqb1f3+/l49KNTXTVVHzrMp9+p/lZ9LYSAiIjUXi52ASIizUoNWEQkEjVgEZFI1IBFRCJRAxYRiUQNWEQkEjVgkSZhZtub2Qux65AiNWBpema2n5lNN7OPzewtM7vbzLbr5z5fM7NdqlVjBccLZrbesrYJITwcQtiwj/t/zcwWmdmILq/PTI7d2pf9Njs1YGlqZnY8MBk4BxgJrAVcCkyMWVe1mVlLFXbzKrBvyT7HAkOqsN+mpQYsTcvMVgR+AhwRQvh9COGTEMLiEMIdIYR/TbYZZGaTzezN5DHZzAYl60aY2Z1m9r6Z/d3MHjaznJldgzfyO5JUfVKZY+9oZnPM7CQzm5ck76+b2e5m9mKyv9NKtv+SmT2WHOstM/svMxuYrHso2ez/kuN9u2T/J5vZ28CVhdeSn/l8coxxyfLqZjbfzHZcxlt2DXBAyfIk4Oo+vfkCqAFLc9saGAzcuoxtTge2AjYHNgO+BPwoWXcCMAdYFU/PpwEhhPBd4HXgayGEoSGEn3ez79WS468BnAH8GvgOsAWwPfBvZrZOsm0bcBwwIql7AnA4fsAdkm02S453Y8n+VwbWBg4uPXAI4RXgZOBaMxsCXAlcFUJ4cBnvxVRgBTP7JzPLA/sA1y5je+mBGrA0s1WAd0MIS5axzf7AT0II80II84GzgO8m6xYDo4C1k+T8cOjdzVUWA2eHEBYDv8Wb60UhhI9CCM8Az+JNnxDCjBDC1BDCkhDCa8CvgH/uYf/twI9DCJ+FEBZ2XRlC+DXwMjAt+e84vYKaCyn4y8BzwNwKfka6oQYszew9YEQP50dXB2aXLM9OXgP4Bd7A/sfMZpnZKb09fgihLflzoUG+U7J+ITAUwMw2SE53vG1mH+LnrDt9IVbG/BDCP3rY5tfAGODiEMJnFdR8DbAf8D10+qHf1IClmT0GfAZ8fRnbvIn/E75greQ1kqR6QghhXeBfgOPNbEKyXbVvM3gZ8DywfghhBfx0h/XwM8uswcyG4l9AXgGcaWYr91RECGE2/mXc7sDvK6hblkENWJpWCOED/NzrJckXYEPMbICZfdXMCudtbwB+ZGarJkOwziA572lme5rZemZmwAf4edr25OfeAdatYrnDgA+Bj81sI+CwLuv7cryLgOkhhB8CfwR+WeHPHQjsHEL4pJfHky7UgKWphRD+Azge/2JtPvAGcCRwW7LJvwPTgSeBp4DHk9cA1gfuAz7G0/SlIYQHknXn4o37fTM7sQqlnoj/0/8j/LTBjV3WnwlclRxv7552ZmYTgd0oNvLjgXFmtn9PPxtCeCWEML0XtUs3TDdkFxGJQwlYRCQSNWARkUjUgEVEIlEDFhGJpBo36JAGMGLEiNDa2hq7DJFMmTFjxrshhFX7+vNqwAJAa2sr06drZJFIb5jZ7J636p5OQYiIRKIGLNJMPvgAHn0UZs8GXQMQnRqwSKMLAW66CcaMgeHDYdttobUVPvc5OOUU+Oij2BU2LTVgkUb29tswYQLsvTfk83D22XD77XDppbDTTvCzn8GGG8K998autCnpSziRRjVnjjffOXPgkkvgkEO8CRccdhhMmwYHHQRf+xrccAPstVe8epuQErBII5o7F7bf3hPw//wPHH545+ZbsOWW8PDD8MUvekq+6aba19rE1IBFGs3ixd5M330X7r/fz/kuy4or+imIrbeGSZPg6adrU6eoAYs0nJNO8pEOV1wB48dX9jNDh8LNN3sz3msv+PDDdGsUQA1YpLHceSdMngxHH+0puDdWWw1uvBFeeQWOOCKd+qQTNWCRRvHxx36ud8wY+MUv+raPHXaA006Da6+FP/2puvXJUtSARRrFmWfCG2/Ar34FAwf2fT+nnQbrr+/N/B89zekp/aEGLNIInnzSTz0cdBBss03/9jV4sI8TfvllOPfc6tQnZakBizSCk06CFVaA886rzv522QW+/W0/lfHmm9XZpyxFDVgk6x54wIeRnXYarNzjzPKVO+ccWLIEfvKT6u1TOlEDFsmyEPx+DqNHw5FHVnff667rV89dfjm8+GJ19y2AGrBItt16K/z1r3DWWX7uttp+9CPf77/9W/X3LWrAIpkVAvz7v/uIhQMOSOcYI0fCscf6JcrPPZfOMZqYGrBIVt19N8ycCaeeCi0p3lfr2GNhueWq9wWfdFADFsmiQvpday34znfSPdaIEXDwwXDddTBrVrrHajJqwCJZ9Oc/w2OP+fCzAQPSP96JJ/rd1H7+8/SP1UTUgEWy6Oc/9xktfvCD2hxvjTX8Tmn//d8wb15tjtkE1IBFsubZZ/3875FH+rnZWjn+ePjsM79KTqpCDVgkayZP9qFhhx5a2+NutBHssYc34IULa3vsBqUGLJIl8+bB1Vf7sLNVV6398U84AebP9y/kpN/UgEWy5Je/9NMAxx0X5/g77gibbw4XXKBp7atADVgkKxYtgssug69+1U8HxGDmzf+553y6I+kXNWCRrLj5Zp9k86ij4tbx7W/76Y+LL45bRwNQAxbJiosv9suOd901bh2DBvmFGXfcAa++GreWjFMDFsmC6dNh6lQfeparg1/bww7zOjQkrV/q4G9SRHp08cU+c/H3vhe7ErfGGj578uWXw6efxq4ms9SARerdu+/6bMUHHOCzXtSLI4+E99+HG26IXUlmqQGL1LsrrvChZ4cfHruSzrbbDsaOhUsu0ZC0PlIDFqlnbW0+9GzHHWGTTWJX05mZfyjMnOnnp6XX1IBF6tldd8Hs2XDEEbErKe873/HTIv/1X7ErySQ1YJF6dsklsPrqMHFi7ErKGzrU75J2003wzjuxq8kcNWCRevXSSz7b8SGH1Oaev311+OGweLGPiJBeUQMWqVeXXeZTDR10UOxKlm2jjWDCBPjVr3wae6mYGrBIPfr0U7jySh9rO2pU7Gp6dsQR8MYbcOedsSvJFDVgkXp0/fU+xrZev3zr6mtfgzXX9HPWUjE1YJF6E4I3srFjfaxtFrS0+A3i77sPnn8+djWZoQYsUm8eeQSeeMKvNDOLXU3lfvhDGDhQQ9J6QQ1YpN5cfDEMHw777x+7kt753Odgn33gqqvgww9jV5MJasAi9WTuXLjlFjjwQFh++djV9N5RR8HHH/vsydIjNWCRevLLX0J7e/3d96FS48fDllv6aYj29tjV1D01YJF68Y9/+FjaPfeEddeNXU3fHXWUX0Ryzz2xK6l7asAi9eK663zG4WOPjV1J/3zrWz52+cILY1dS99SARepBCN6wNt0UdtopdjX9M3Cgj+C47z546qnY1dQ1NWCRenDfffDMMz7jcJaGnnXnkENgueXgootiV1LX1IBF6sGFF8LIkbDvvrErqY5VVvEZPK69FubNi11N3VIDFontqafg7rv9suNBg2JXUz3HHQeLFmn6+mVQAxaJ7ec/9zG/WbnvQ6U23NDvY3zJJT42WJaiBiwS0+zZPqnlQQfByivHrqb6Tj4ZFiyAX/86diV1SQ1YJKYLLvAv3Y4/PnYl6dhqK/jnf/b/zkWLYldTd9SARWKZN8+T4f77+60cG9XJJ8OcOXDNNbErqTtqwCKx/OIXPt38qafGriRdu+0GW2wBZ5/tUxdJBzVgkRjmzYNLL4X99vMvqxqZGZx5Jrz6qlJwF2rAIjGcf77f++FHP4pdSW3ssYdScBlqwCK19vbbPjRr330bP/0WFFLwrFm6VWUJNWCRWjvrLB8RcOaZsSuprT328FERZ57pk46KGrBITb34oo98OOQQWG+92NXUlplfdPLmm7pHREINWKSWTjvNb1JzxhmxK4lj++19BuXzzoP33otdTXRqwCK18tBDPt3QiSf6/GnN6txz/dLkH/84diXRqQGL1MKSJX6P3LXXhn/919jVxLXJJj7l0mWX+ezPTUwNWKQWLrnE73p24YUwZEjsauL76U/9lpVHHNHUc8epAYuk7c03/ZzvV74CX/967Grqw/Dh8LOfwaOP+jT2TUoNWCRNIfiIh0WLfKbgRpjtolomTfIv5Y47DubOjV1NFGrAImm69lq480445xxYf/3Y1dSXXA5+8xv/cDr4YP+wajJqwCJpmTMHjj4att3Wn2Vp663noyLuusubcZNRAxZJw+LFsM8+Pvrhyishn49dUf066iifCfqoo3xi0iaiBiyShtNPh0ce8avedOph2XI5uO46WGEF+Na3mmr6IjVgkWq76Sa/1++hh3oKlp6NGgXXXw/PPw8/+EHTDE1TAxappkcfhe9+F7bZxsf8SuV23tnvFXHTTf4viCbQErsAkYbxwgs+C/Caa8Ltt8PgwbEryp4TToBXXvF7Ray1Fhx2WOyKUqUGLFINzz/vCc7Mv9EfMSJ2RdlkBhdf7CNIDj/cv7w8+ODYVaVGpyBE+uvpp/1b/PZ2ePBBfenWXy0tcPPNfv/gQw7xC1galBqwSH/cfbef7zWDBx6AjTeOXVFjGDTI7xw3caIPTzvuOGhri11V1akBi/RFW5tf3bbnnn4xwV//Cv/0T7GraiyFJnzssTB5Muy+u0/n1EDUgEV6a9YsP+Vw+uk+bvWhh2D06NhVNaZ83keTTJni7/PYsXDbbbGrqho1YJFKffKJ39Vs4439PrZXXw033ABDh8aurPEddBDMmOEfdP/v//n54RdfjF1Vv6kBi/Tkgw98WNQ66/h9bL/xDXj2WR/vq7ub1c7GG8O0aXD++fDww758wAH+d5FRasAi5SxZAv/7v/C97/lVWqeeClts4ZcXX3+9TjnEMnCgjxV+8UU45hg/R7zJJrDDDj7d/YIFsSvsFQtNeAs4Wdr48ePD9OnTY5cRT1ubX0jxyCM+muHee+Hvf/fTC/vt58Ohxo2LXaV0NX++30XtiivgpZf8nPEOO8CECf48bhwsv3xqhzezGSGE8X3+eTVggSZowIsWeTqaPx/eeccH+s+eDS+/7I336afh009929VWg1128XONu+6a6i+wVEkI8Le/wa23wh//6NM/gZ8i2mADH6Gy/vrQ2upXKo4a5ROjrryy//328VSSGrBUxfhVVw3TJ05M/0Dd/f9WeL10fQhLP9rb/bmtrfhYssQfixfDZ5/5Y+FCf3z8MXz0kf+5nNGjYcMNYcwYT0tbbum/sDq3m23vvef35Zg5078wfeEF/7BdtGjpbVtaYNgwfwwZAsst50PgBg3yUx4DBvg2+XzxsfrqcMEFasBSHeMHDgzTazVVenfNrfB66Xqzzo9czp9Lfxnyef8lGTCg+Iuz3HL+GDbMTyMMH+6PVVeFkSO98Y4e7dtKc2hv93/9vP66P8+b5/8qWrDAP6Q/+sj/FbRwYfGDfNEi/2BfsqT4gd/e7rNb33tvvxuw7gUhbtNNoZFPQYjkcn7qYdSo2JV00CgIEZFI1IBFRCLROWABwMw+Al6IXUcPRgDvxi6iB1moEbJRZxZq3DCEMKyvP6xzwFLwQn++TKgFM5uuGqsjC3Vmpcb+/LxOQYiIRKIGLCISiRqwFEyJXUAFVGP1ZKHOhq9RX8KJiESiBCwiEokasIhIJGrATc7MdjOzF8zsZTM7JXY9BWa2ppk9YGbPmtkzZnZM8vrKZvYnM3speV6pDmrNm9lMM7szWV7HzKYl7+mNZjYwcn3DzexmM3vezJ4zs63r7X00s+OSv+enzewGMxtcD++jmf3GzOaZ2dMlr5V978z9Z1Lvk2bW4/1L1YCbmJnlgUuArwIbA/uaWb1M67sEOCGEsDGwFXBEUtspwP0hhPWB+5Pl2I4BnitZ/hlwYQhhPWABcGCUqoouAu4JIWwEbIbXWjfvo5mtARwNjA8hjAHywD7Ux/v438BuXV7r7r37KrB+8jgYuKzHvYcQ9GjSB7A1cG/J8qnAqbHr6qbW24Ev41frjUpeG4VfQBKzrtHJL+HOwJ2A4VdvtZR7jyPUtyLwKskX7iWv1837CKwBvAGsjF8cdiewa728j0Ar8HRP7x3wK2Dfctt191ACbm6F//EL5iSv1RUzawW+AEwDRoYQ3kpWvQ2MjFRWwWTgJKA9WV4FeD+EsCRZjv2ergPMB65MTpNcbmbLU0fvYwhhLnA+8DrwFvABMIP6eh9Ldffe9fr3SQ1Y6pqZDQVuAY4NIXxYui54zIg2jtLM9gTmhRBmxKqhAi3AOOCyEMIXgE/ocrqhDt7HlYCJ+IfF6sDyLP3P/rrU3/dODbi5zQXWLFkenbxWF8xsAN58rwsh/D55+R0zG5WsHwXMi1UfsC3wL2b2GvBb/DTERcBwMyvcZyX2ezoHmBNCmJYs34w35Hp6H3cBXg0hzA8hLAZ+j7+39fQ+luruvev175MacHP7G7B+8m3zQPyLjz9Ergnwb5SBK4DnQggXlKz6AzAp+fMk/NxwFCGEU0MIo0MIrfh7978hhP2BB4BvJpvFrvFt4A0z2zB5aQLwLHX0PuKnHrYysyHJ33uhxrp5H7vo7r37A3BAMhpiK+CDklMV5cU68a5HfTyA3YEXgVeA02PXU1LXdvg/7Z4Enkgeu+PnWO8HXgLuA1aOXWtS747Ancmf1wX+CrwM3AQMilzb5sD05L28DVip3t5H4CzgeeBp4BpgUD28j8AN+Hnpxfi/Jg7s7r3Dv4C9JPldegof1bHM/ffrUmQz2w3/J1ceuDyEcF6fdyYi0mT63ICTMaQv4kOD5uD/nN03hPBs9coTEWlc/bkh+5eAl0MIswDM7Lf4N5ndNuARI0aE1tbWfhxS+uv9933y1zXX7Pz6jBkz3g0hrFpY/nLuW73/ZO4627Hluix2s77L61bYTy7Xeb/JcnF9YRblLvvJ5TuWO7bN5zvvq+P1ws/6c8h1OXa+cw2h4/XulpPn5OdC3sqvT57bW7r+XOGZzsvJYdq7rO+6XNiuuL7zftpbOq9f6rlwnK7btYRu9hs6r0+eSV6nJWDJn63FR+rl88lzsjxggI80G5BvA2Bgiz8Pyheeff1yLYsBGJw8L5/3KeaXyyfLLZ8BMCTnrw/L/wOAocnzCrmFyfrPkmV/fVjHs+9nmIVk2d+EobnBAORWe6mb6bz7rj8NuNyYty27bmRmB+NXhbDWWmsxXTPvRnXSSXDxxUtPgGxms+NUJNK8Up+SKIQwheSemePHj9e9LyNbvBgGDEhp54XTWYV0GZJrE5KEGtqTJJTrsr69c4ItnBaz9mR9IWUmy4XUaYVLH3Jd9kNb8pzvSHTWlrxWSMIFbe2dFi0ZGBTo/HohCRdqKlzDZHRdLugoLllPl/XJ2uQyg/alfhMLW3b+yVyy3N7NcleFd6Q92S6XbNdedusydXVTT3G/5Y8byv7Zf6qNrqrchvq7uyQJ0164BiRJ0P3cbdlD9eNn63oMqZS3aBEMjHprGBEp6M9nRccYUrzx7gPsV5WqJDWffgrLL5/yQeoqCeeTbZMSlISTZSXhHnVJwmkk4D6XGEJYYmZHAvfi/5f/JoTwTNUqk1R8+CEMHRq7ChGBfn5GhBDuAu6qUi1SA++9B6usUqOD1UUSLp4P9m2TEpSEk+X0k3DXL34ym4RToEuRm8zcubDaarGrEBGowSgIqR9tbfDaa/DNb/a4aXXFTMJlRkb4tkkJSsLJspJwDErATeTpp2HJEhgzJnYlIgJ1+Zkgafnzn/15hx0iFRAjCS9jjLBvm5SgJJwsKwnXkhJwE7n2Wth446UvQxaROOros0DSNHUq/O1vfhlydDVMwpVcLefbJiUoCSfL1UvClYwR7rzcPElYCbgJLFoEhx0GI0fCAQfErkZECurgM0DSFAKcfjo88QTcdhusUOnlPGbFpJqWGiTh3tw3ApSE00jCvblarvNy4ydhJeAGFoLf/ez88z0BT5wYuyIRKaUE3KA++ACOPhquvhqOOAL+8z/7sJOOZJrhJNyHO6j59kkJSsLJcn+ScO/vG9F5uXGTsBJwA7rtNh/tcO21cMYZ/sVbTn/TInVHCbiBPPYYnHsu3HEHbLqpN+IvfrGPO7NcSRLNcBLux72EffukBCXhZLkvSbjvd1DrvNx4SVi5KOPa2uDWW2HbbWGbbeAvf4FzzvEZL/rcfEWkJpSAM+q55+DGG+G66+Dll6G11c/zfv/7VbzdZGGutQwn4WrMquHbJyUoCSfLlSfhatxLuPNy4yRhNeAMefFFb7q/+53f18HMLys++2z4xjegRX+bIpmiX9k69vHH8PDDcN998Kc/wVNPedPdbjv/Ym2vvWDUqHSObTnrSJpZTsLVnF/Ot09KUBJOlntOwtWcVaPzcvaTsBpwHVm8GKZN84Z7//1++fCSJT6H27bbwuTJfivJNdaIXamIVIMacETz5nnDLTweeww++cQD3xZbwAknwC67ePNdbrna11dIlplOwinMtOzbJyUoCSfLTZCEU1D/FTaIzz7zy4GnTvVmO3UqvPqqr8vnYexYmDQJJkyAHXeElVeOWq6I1IAacAo+/hiefNIb7syZ/vzkk35THPBTCFtt5ZcHb7mlp93UZyrurZJxwJlOwqnMtAxKwqVHVxLuq/qtLCPeeadzo505E156qdgrVl4ZvvAFOOYYb7ZbbgmjR8etWUTqgxpwhRYuhGef9ZEIpY+33y5u09rqzXb//WHzzf3Po0cXA1um5KyY+rKchFOZaRmUhLvsv6MaJeHeqL+KImtrg1mzlm60L7/c8bvJ4MF+r4Vddy022s02g+HD49YuItnS1A143rylG+0zz8Cnn/p6M/j85/0Lsn328eexY2G99ZYOOA2pcC41w0k4lZmWS/ejJNx5/x3VFJNwf+4l3HmfjZeE66eSFH36qTfWrs123rziNquu6s31oIOKjXaTTerwyzERaRgN14D//nf/Iuzxx4vPL75YDEuDB/u07HvsUWy0Y8f6dD1SZFa8iizLSTiVmZZBSbgXSbgWMy2XVrz0cv0m4fgV9FEI8OabSzfb118vbrPmmjBuXOfTB5//fJOcPhCRupeZBtze7qcRHnrI74/w0EPw1lvF9RtsAFtv7bM/fOEL/hgxIl69mZfLFVNXlpNwCjMt+/rkmErCnZa7ypHOTMvLOmaWknDdNuDFi2HGjGKzfeQRWLDA162xhl8tttVWnnA32wyGDYtarohIr9VVAw7Bbyh+9dV+y8UPP/TXN9jAb7e4/fZ++8XW1oyOrc0Ss470l+kknMJMy6Ak3Jsk3Jcxwp1q6qaWRkjCPR7RzNYErgZG4v8NU0IIF5nZysCNQCvwGrB3CGFBX4qYNQuuucYb76xZPvJgr73ga1/zWy+utlpf9ioiUt8qaflLgBNCCI+b2TBghpn9CfgecH8I4TwzOwU4BTi5twX8x3/AiSd6WJkwAc4809Ouhn9FVjKjcKaTcBozLYOSsJJwVfR4pBDCW8BbyZ8/MrPngDWAicCOyWZXAQ/SywY8ZYo33732ggsv9FELIiLNolet3sxagS8A04CRSXMGeBs/RVGxV1+FQw+F8ePh+uv9puNSP8ysI91lOgmnMdMyKAkrCVdFrudNnJkNBW4Bjg0hfFi6Lvj/nWV/i8zsYDObbmbT58+f3/H66NGw004+fveuu/pWvIhIllXU4s1sAN58rwsh/D55+R0zGxVCeMvMRgHzyv1sCGEKMAVg/PjxHU16wAC47Tb48pf9FMSuu/oNySdO9KvVJLKcdaQ5JWGUhPuRhNOYabnzdtlNwj0mYDMz4ArguRDCBSWr/gBMSv48Cbi9twcfNgzuuQdOPtnvzbDPPj7i4eCDfdxve09/EyIiGVZJa98W+C7wlJk9kbx2GnAe8DszOxCYDezdlwKGD4dzzoGf/hQefBCuugquuw5+/WtYaSUfhrbDDj4GeNw4T85SA5brSG9KwigJ9yMJp3Ev4U41dVNL9ZNw9VUyCuIvLP33UTChWoXk8z4MbcIEuPRSuP12eOABvxLujjt8myFD/HLjwgUZX/qShquJSHbV1ZVwBUOH+qwS++/vy2+/7VfIFe4DcdZZxSGeG27o930YN654DwhNaFkFpTNiKAkrCRfW9iEJpzmrRqeauqmlnpNwXTbgrlZbDb75TX8AvP8+PPoo/O1vPoriL3+BG24obr/22sWGPG6cz1qx+uq6fFlE6ksmGnBXw4fD7rv7o+Ddd70ZFx6PP+6jLApBZ6WVirek3HRTfx4zRjfx6VYuT8dnv5KwknDnvfYqCddifrlONXVTS3+TcBoy2YDLGTHCh7R9+cvF1z76CP7v//zx1FM+NfzVV/vrBa2tnW/MPnas3/yn6b/syxmFZqBGrEbcv0bc+4s1ikeup0ZcfQ3TgMsZNsxHUWy3XfG1EGD27KWnJ7r7bliS/MUNHAgbbbR0Y87sDMciUpcaugGXY+apt7XV77ZW8Nln8MILnpILTfnPf/YhcQXDh/tpi66nMlZYodb/FenzS5ELS0rCSsL9ScJ9v2y5eOT4STgNTdeAuzNokDfUTTft/PqCBfD0053T8nXXFe9VDD5LcmF6+sLzaqspLYvIsqkB92CllXzc8fbbF18LAd54w9PyE0/44/HH4eabi9t87nPFhrz55rDFFt6oM9OU8/mOtKUkjJJwv5Jw/2/gUzxyYyVhNeA+MIO11vLHnnsWX//gA2/KM2d6U545Ey64wKdXAh+fvOWW/thqK7+QZKWV4vw3iEh8asBVtOKKS6flRYvg2Wdh+nSYNg2mTvX7XxQC1gYbeDPeckvYZhs/BZLLld9/TZl1pCslYZSEO36geZNwGtSAUzZwYPE0xA9/6K99+KE35KlTvSnfc48PjwMfTrfTTrDLLn5Z9rrrZui0hYj0ihpwBCusADvv7A8oDo176CG4/35/3HSTr2ttLd4jY7fdanjKorTrKwkrCffPc1UpAAAO+klEQVQjCac55b2vz24SVgOuA6VD4w44wH/XX3ih2IxvuQWuuMIvDvnKV2Dvvf2+ySuuGLtyEekPNeA6ZOYXgmy0ERxxBLS1+X0vbrkFfvc7+OMf/dTGbrt5M95rrxRuYp/PLZWelIRREu5DEk5zeqPikbOZhOvh6x7pQT7vX9T94hfw2mvw2GPemGfMgO98x28+9NOfwnvvxa5URHpDDThjzLwZX3ABvP463HefjzE+4wyfVfrII/18cr/lcp6G8jk/aOkjn/dxwmaYmae6nPkNfHL54rLl/JEsF7fP+aOwv2S5Y31HDV32U3gPclZMmP5Cp/Ud+01TCJ2Tdmgvpl08CXck9dL17cEfHbsJnobb2/1R2G+yXFyfPLrup70tefhyx/Ztbf4o7K/waGv3R7J/aw9Ye8nxC+uT7a293dNwW4C2csvJo609eQQsWddpfXsgt8QfxZ8pPEgeyXK7J/9cWyBXsr7rcmG74np/FPaTW+Ipt7j/Lo/Ccbput8T80WW/aVADzrBczr+cu+uu4pROU6bAxhvD5Mn++yci9UsNuEGMGQO/+Q289BLsuCMcd5yPK37mmb7tL+RKkmqWk3DalIR7TMKlKTjLSTgNasANZu214c474frrYdYsb8KPPBK7KhEpR6MgGpAZ7Luv34Zzl1186Nof/uCnKyqWy3V8Q25dP6ezNDqiFiMjSvev0RGJ4uiIvt1Bzfdauhx/dET1KQE3sDXX9Is71l3Xp3OaOzd2RSJSSgm4wY0cCbfe6veYOOggH0Nc0anRfPGzOdNJuJZjhEv3ryScaKdv940o3bJeknD1KQE3gfXWg3PO8Vk/HnggdjUiUqAG3CQOPdRvh3nppRX+gFlx9EMyOiHkLHujI2KMEQaNjigdAVHBGOEsjI5Igxpwkxg8GCZNgttvh08+iV2NiIAacFP5yld84tGpU3veNpQm0Awn4ZIFJeFISbg3V8vVcxJOgxpwE9l6a3/+61/j1iEiTqMgmsiKK/qoiFdeqWDjvBW/0S58E5/v/HmdhdERdXEHtdL9N+HoiHRmWl76J2s1OqKalICbTGur38RHROJTAm4yI0bAW2/1vF3I5ZYe25nFJFxP9xIu3X8TJeF05pcr3TK7SVgJuMmstBIsWBC7ChGBXiRgM8sD04G5IYQ9zWwd4LfAKsAM4LshhEXplCnVMmQILFxYwYal54CznITrcVaN0v03QRJOd6bl0i2zl4R7k4CPAZ4rWf4ZcGEIYT1gAXBgNQuTdAweXGEDFpHUVZSAzWw0sAdwNnC8+SDMnYH9kk2uAs4ELkuhRqmiAQNg8eKetws5K8kVGU7C9Ty/XOn+GzgJpzHTMtQ+Caeh0gQ8GTiJ4ju0CvB+CKFwgd4cYI1yP2hmB5vZdDObPn/+/H4VK/3X0uIXY4hIfD0mYDPbE5gXQphhZjv29gAhhCnAFIDx48enHCekJ7lcMTAuS8jnoCPJJq9lMAlnYqbl0v03YhJOYaZlqH0STkMlpyC2Bf7FzHYHBgMrABcBw82sJUnBowHdbTYDKm3AIpK+HhtwCOFU4FSAJAGfGELY38xuAr6Jj4SYBNyeYp1SJWaVhbmQN0ozAGQzCVdlVg1QElYSTkV/9n0y/oXcy/g54SuqU5KkqdIGLCLp69WVcCGEB4EHkz/PAr5U/ZIkTRXfACxvJWe+MpyEqzm/HCgJ9yEJV3N+OV9Pl/XJ2pSTcBp0JZyISCS6F0STqTQBdx4HXJDBJJzGTMugJNyLJJzGTMu+ni7rk7UZSsJKwCIikSgBS1mhZDaJLCfhVGZaLlmvJJyUoSTcJ0rAIiKRKAFLWe0tttRMsFlMwtWYVQOUhPuVhKswq0bn5YLaJuE0KAGLiESiBCxlhZx1JIZMJ+Eqzi8HSsLNnITToAQsIhKJErCU5feCcJlOwinMtAxKwr1KwinMtNx5uSDdJJwGJWARkUiUgKWskIeu2SGbSTiFmZZBSbgXSTiNmZZ9+9om4TQoAYuIRKIELGX5OeDyCSBbSTiFmZZL96MknKzvPgmnMdOyr691Eq4+JWARkUiUgKUsTy3LPheWhSScykzLoCTcqyRc/ZmWS2uqXRKuPiVgEZFIlIClrPa8lYx/zG4STmOmZVAS7l0STmF+OYiQhKtPCVhEJBIlYCkr5MvdDSqDSTiFmZZLa1ISriAJpznTMmQ6CSsBi4hEogQsZZWeA1YSVhLuVxJOYaZlqH0SToMSsIhIJErAUla5c8CZTMJpzLQMSsK9SMJpzC8HtU/CaVACFhGJRAlYygolH81KwigJ9yMJpznTMmQ7CSsBi4hEogQsZYX80q9lMQmnMdOyHzM5hpKwPy8rCacxvxzUPAmnQQlYRCQSJWApqz3f/adzlpJwKjMtg5Jwb5JwmjMtQ82ScBqUgEVEIqkoAZvZcOByYAweUH4AvADcCLQCrwF7hxAWpFKl1FzIG+1JFlUSBiXhvifhVGZahghJuPoqTcAXAfeEEDYCNgOeA04B7g8hrA/cnyyLiEiFekzAZrYisAPwPYAQwiJgkZlNBHZMNrsKeBA4OY0ipfbaWyCXZIBMJ+FUZloGJeHKk3AaMy1DjCRcfZUk4HWA+cCVZjbTzC43s+WBkSGEt5Jt3gZGlvthMzvYzKab2fT58+dXp2oRkQZQyTngFmAccFQIYZqZXUSX0w0hhGBmZT+KQwhTgCkA48ePT/njWqrF7wXhspyE05lpGZSEe5GEU5hp2bdPSqhREk5DJXueA8wJIUxLlm/GG/I7ZjYKIHmel06JIiKNqccEHEJ428zeMLMNQwgvABOAZ5PHJOC85Pn2VCuVmiq9Ei7LSTiNmZZBSbhXSTiVmZah1kk4DZVeiHEUcJ2ZDQRmAd/H/+//nZkdCMwG9k6nRBGRxlRRAw4hPAGML7NqQnXLkXpR/l4QLktJOI2ZlkFJuFdJOJWZlqHWSTgNuhJORCQS3QtCygrL+GjOUhJOZX45UBIuqCAJpzLTcul+apSE06AELCISiRKwlNXe0vNssFlIwlWdVQOUhLuzjCScykzLUPMknAYlYBGRSJSApaxOV8JlOgn34r4RoCTcX+WScAozLfv65Ji1SsIpUAIWEYlECVjKCi0BOhKty2YSrvy+EaAkXDWlSTiFmZahMZKwErCISCRKwFKWXwnXOXlmMQmnMdOyb68kXJEQenffCKjbJJwGJWARkUiUgKWszrMiZzcJV2VWDVAS7o80ZlqGmifhNCgBi4hEogQsZYV8KEmgHa8mz0rCdGyvJFwxJeGlKAGLiESiBCxl+Thgl+UknMpMy6Ak3B9ZTcIpUAIWEYlECVjKKk3ABVlMwmnMtAxKwlWhJKwELCISixKwlJcPdJd9lISVhKsqK0k4BUrAIiKRKAFLeSXngLOdhPtzL+FidUrCSsJpUAOWsqzMKYhsNuLqTW/Ueb0acWrqtBGnQacgREQiUQKWsqylncLnc7aTcN9vZVn6U0rCSsJpUAIWEYlECVjKyufbSz73s5uE05zyvvN6JeHU1E0Srj4lYBGRSJSApax8SzGJZTkJpzvlfbE6JeEmSMIpUAIWEYmkogRsZscBP8QDwFPA94FRwG+BVYAZwHdDCItSqlNqbMCAJXT93yOLSTjdKe9Lj64k7MsNnIRT0GMCNrM1gKOB8SGEMfjf9D7Az4ALQwjrAQuAA9MrU0Sk8VR6DrgFWM7MFgNDgLeAnYH9kvVXAWcCl1W7QIljQL507GN2k3CaU96X/pSScOMn4TT0uOcQwlzgfOB1vPF+gJ9yeD+EUPjfcg6wRrmfN7ODzWy6mU2fP39+daoWEWkAPSZgM1sJmAisA7wP3ATsVukBQghTgCkA48ePT/HjUKppYEu5q3+UhJWE+5CE00zBXmTxWJ2OXd0knIZKsvUuwKshhPkhhMXA74FtgeFmVviNHA3MTalGEZGGVMk54NeBrcxsCLAQmABMBx4AvomPhJgE3J5WkVJ7g/LLuv5dSVhJuBdJuBbng0v3n1YSTkEl54CnATcDj+ND0HL4KYWTgePN7GV8KNoVqVUpItKAKhoFEUL4MfDjLi/PAr5U9YqkLgzK9xA7gSwk4d7cN6J0WUm4ikm4liMjSvdf7SScAl0JJyISie4FIWUt17K4F1vXbxJOc6JPUBKuKAnHGCNcuv8qJeE0KAGLiESiBCxlDe5VAi6ovyRciynvQUl4WUk46tVypfvvZxJOgxKwiEgkSsBS1vL5/tzYrn6ScF/uoObrlYT9uMlx+pGE6+K+EaX772sSToESsIhIJErAUtZy+b6cA+5KSVhJuDgOOPNJOAVKwCIikSgBS1nLt3xWxb0pCTd1Eq7HewmX7r/SJJwCJWARkUiUgKWsIblFKfzfoSTclEm4nmfVKN1/D0k4DUrAIiKRKAFLWcPy/yguZDgJpzHTsq9XEvbjJsdZRhLOxPxypfvvJgmnQQlYRCQSJWApa2hpAi7IYBLuz72EO++zu2MqCftxk+OUScKZmmm5dP9dk3AKlIBFRCJRApayVsgt7H5lppJw/2fV6LzP7o6pJOzHTY5TmoTTmGnZd0SquibhFCgBi4hEogQsZQ3JVXAlXCaScPXmlyvuc1nHVBL24ybHyeXSmWkZap+EU6AELCISiRKwlLVCrswoiO7UcRLuz9VynbZTEvbt+pCEU5lpuWR9zZJwCpSARUQiUQKWsob1JgEX1GESTmOmZV9WEobKknAqMy37C53WZzEJKwGLiESiBCxlDcv1Y0YMJeEyx23iJJzCTMu+mP0krAQsIhKJErCUNcwC9CcFg5Jw2eM2YRJOYaZlaIwkrAQsIhKJErCUNSzXAu1J3MpwEk5zVo1O2ykJ+3Zlk3AKMy1DQyRhJWARkUiUgKWsobnBQDIWWElYSbhfSTiFmZZL95PhJKwE3GS22gqOOSZ2FSICYKGGnwZmNh+YXbMDSm+sHUJYNXYRIs2kpg1YRESKdApCRCQSNWARkUjUgEVEIlEDFhGJRA1YRCQSNWARkUjUgEVEIlEDFhGJRA1YRCQSNWARkUjUgEVEIlEDFhGJRA1YRCQSNWARkUjUgEVEIlEDFhGJRA1YRCQSNWARkUjUgEVEIlEDFhGJRA1YRCQSNWARkUj+P4VunBEI1yzqAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 360x360 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#%% plot the distributions\n",
    "\n",
    "pl.figure(1, figsize=(6.4, 3))\n",
    "pl.plot(x, a, 'b', label='Source distribution')\n",
    "pl.plot(x, b, 'r', label='Target distribution')\n",
    "pl.legend()\n",
    "\n",
    "#%% plot distributions and loss matrix\n",
    "\n",
    "pl.figure(2, figsize=(5, 5))\n",
    "ot.plot.plot1D_mat(a, b, M, 'Cost matrix M')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solve EMD\n",
    "---------\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFgCAYAAACFYaNMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl4VdW5x/Hvm4RBJgFBQECj4nAVUDEtiENVrFqHi7dYxeoVby04D6h1bK1tr9apiperVipaB7SKVmmpw1XUakXQUBwYZBBEQTGACDIIJFn3j3eniRDIdJJ1ht/nefZzcvbZ55w3B/LLytprr2UhBEREpOnlxS5ARCRXKYBFRCJRAIuIRKIAFhGJRAEsIhKJAlhEJBIFsEiWMLNDzWxO7Dqk9hTAIgkzO8vMPjCzdWa21MzuNbP2yWO/N7M1ybbRzDZVuf98E9QWzKzXto4JIbwRQtirAe8x1MymmtlaMytJvj7fzCx53MzsFjNbkWy3VDwm9aMAFgHM7HLgFuBnwPbAAGAX4CUzax5CODeE0CaE0Aa4CXii4n4I4QfxKndmVtDA518O3AXcBnQFugDnAgcDzZPDRgAnAfsBfYETgXMa8r65TgEsOc/M2gG/Ai4KIbwQQtgUQvgYOAUoBM6ox2sebmaLzezKpDX5uZmdZGbHmdlcM/vSzK6tcvx3zewtM/sqOfZ/zax58tjryWHvJS3uU6u8/lVmthR4sGJf8pzdk/fol9zfycyWmdnh1dS6PfBr4PwQwlMhhK+Dmx5COD2EsCE5dBjwuxDC4hDCEuB3wFl1/WykkgJYBAYCLYE/V90ZQlgDPAd8v56v2zV53e7A9cAf8DA/EDgU+IWZ7ZocWwaMBDoBBwGDgPOTOg5LjtkvaXE/UeX1O+It9RGb1f4RcBXwqJm1Ah4EHgohvFZNnQcBLYAJNXw/+wLvVbn/XrJP6kkBLOKhtzyEUFrNY58nj9fHJuDGEMIm4E/J69yVtDBnArPwP+cJIUwLIUwJIZQmre/7gO/V8PrlwC9DCBtCCOs3fzCE8AdgPjAV6AZct5XX2eL7N7PJSWt8vZlV/AJoA6yq8rxVQBv1A9efAlgElgOdttKP2i15vD5WhBDKkq8rAvKLKo+vx0MNM9vTzCYmJ/9W4/3MNQX/shDCNzUc8wegNzC6SlfCFnWy2fcfQhgYQmifPFaRE2uAdlWe1w5YEzSjV70pgEXgLWAD8MOqO82sDfADYFIT1HAv8CGwRwihHXAtUFPLcpvBl9Q/ChgL3GBmHbdyaMX3P7iG95tJ0mJP7Jfsk3pSAEvOCyGswk/CjTazY82smZkVAk8Ci4FHmqCMtsBqYI2Z7Q2ct9njXwC71fE17wKKQwg/Bf4G/L66g0IIX+Hf/z1mdrKZtTWzPDPbH2hd5dCHgcvMrLuZ7QRcDvyxjjVJFQ0auiKSLUIIt5rZCuB2YHc8DJ8FTt/Gn+6pdAUwBrgSmA48ARxZ5fEbgIfMbDv8hFvJtl7MzAYDxwJ9kl2XAe+a2ekhhHGbH598/0uS938YWAsswE/kTU4Ouw//JfBBcv/+ZJ/Uk6n7RkQkDnVBiIhEogAWEYlEASwiEokCWEQkEo2CEAA6deoUCgsLY5chklGmTZu2PITQub7PVwALAIWFhRQXF8cuQySjmNmihjxfXRAiIpEogEVyyapVMHkyLFoEugYgOgWwSLYLAcaPh969oX17OPhgKCyEHXeEq6+Gr7+OXWHOUgCLZLOlS2HQIDjlFMjPhxtvhAkT4J574Igj4JZbYK+94MUXY1eak3QSTiRbLV7s4bt4Mdx9N5xzjodwhfPOg6lTYfhwOPFEePxxGDIkXr05SC1gkWy0ZAkceqi3gP/v/+D8878dvhX694c33oDvfMdbyePHN32tOUwBLJJtNm3yMF2+HCZN8j7fbdl+e++COOggGDYMZsxomjpFASySda680kc6jB0LRUW1e06bNvDUUx7GQ4bA6tWNW6MACmCR7DJxIowaBRdf7K3guujaFZ54Aj76CC64oHHqk29RAItkizVrvK+3d2+47bb6vcZhh8G118Kjj8JLL6W2PtmCAlgkW9xwA3z6Kdx3HzRvXv/XufZa2GMPD/NvalrzUxpCASySDd5/37sehg+HgQMb9lotW/o44fnz4be/TU19Ui0FsEg2uPJKaNcObr45Na931FFw6qnelfHZZ6l5TdmCAlgk0736qg8ju/Za6Li1lefr4aaboLQUfv3r1L2mfIsCWCSTheDzOfToARdemNrX3m03v3ru/vth7tzUvrYACmCRzPbMM/D22/CrX3nfbar9/Of+ur/4RepfWxTAIhkrBPjv//YRC2ee2Tjv0aULXHqpX6I8e3bjvEcOUwCLZKrnn4fp0+Gaa6CgEefVuvRS2G671J3gk39RAItkoorW7847wxlnNO57deoEI0bAuHGwYEHjvleOUQCLZKK//x3eesuHnzVr1vjvd8UVPpvarbc2/nvlEAWwSCa69VZf0eInP2ma9+ve3WdK++MfoaSkad4zByiARTLNrFne/3vhhd4321Quuww2bPCr5CQlFMAimWbUKB8adu65Tfu+e+8Nxx/vAbx+fdO+d5ZSAItkkpISePhhH3bWuXPTv//ll8OyZX5CThpMASySSX7/e+8GGDkyzvsffjjsvz/ccYeWtU8BBbBIpti4Ee69F37wA+8OiMHMw3/2bF/uSBpEASySKZ56yhfZvOiiuHWceqp3f4weHbeOLKAAFskUo0f7ZcfHHBO3jhYt/MKMv/4VFi6MW0uGUwCLZILiYpgyxYee5aXBj+1553kdGpLWIGnwLykiNRo92lcuPuus2JW47t199eT774d162JXk7EUwCLpbvlyX634zDN91Yt0ceGF8NVX8PjjsSvJWApgkXQ3dqwPPTv//NiVfNshh0CfPnD33RqSVk8KYJF0VlbmQ88OPxz23Td2Nd9m5r8Upk/3/mmpMwWwSDp77jlYtAguuCB2JdU74wzvFvnf/41dSUZSAIuks7vvhp12gsGDY1dSvTZtfJa08ePhiy9iV5NxFMAi6WrePF/t+JxzmmbO3/o6/3zYtMlHREidKIBF0tW99/pSQ8OHx65k2/beGwYNgvvu82XspdYUwCLpaN06ePBBH2vbrVvsamp2wQXw6acwcWLsSjKKAlgkHT32mI+xTdeTb5s78UTo2dP7rKXWFMAi6SYED7I+fXysbSYoKPAJ4l9+GT78MHY1GUMBLJJu3nwT3n3XrzQzi11N7f30p9C8uYak1YECWCTdjB4N7dvD6afHrqRudtwRhg6Fhx6C1atjV5MRFMAi6WTJEnj6aTj7bGjdOnY1dXfRRbBmja+eLDVSAIukk9//HsrL02/eh9oqKoL+/b0borw8djVpTwEski6++cbH0p5wAuy2W+xq6u+ii/wikhdeiF1J2lMAi6SLceN8xeFLL41dScP86Ec+dvnOO2NXkvYUwCLpIAQPrL594YgjYlfTMM2b+wiOl1+GDz6IXU1aUwCLpIOXX4aZM33F4UwaerY155wD220Hd90Vu5K0pgAWSQd33gldusBpp8WuJDV22MFX8Hj0USgpiV1N2lIAi8T2wQfw/PN+2XGLFrGrSZ2RI2HjRi1fvw0KYJHYbr3Vx/xmyrwPtbXXXj6P8d13+9hg2YICWCSmRYt8Ucvhw6Fjx9jVpN5VV8HKlfCHP8SuJC0pgEViuuMOP+l22WWxK2kcAwbA977n3+fGjbGrSTsKYJFYSkq8ZXj66T6VY7a66ipYvBgeeSR2JWlHASwSy223+XLz11wTu5LGdeyxcOCBcOONvnSR/IsCWCSGkhK45x748Y/9ZFU2M4MbboCFC9UK3owCWCSG22/3uR9+/vPYlTSN449XK7gaCmCRprZ0qQ/NOu207G/9VqhoBS9YoKkqq1AAizS1X/3KRwTccEPsSprW8cf7qIgbbvBFR0UBLNKk5s71kQ/nnAO9esWupmmZ+UUnn32mOSISCmCRpnTttT5JzfXXx64kjkMP9RWUb74ZVqyIXU10CmCRpvL6677c0BVX+Pppueq3v/VLk3/5y9iVRKcAFmkKpaU+R+4uu8DPfha7mrj23deXXLr3Xl/9OYcpgEWawt13+6xnd94JrVrFria+3/zGp6y84IKcXjtOASzS2D77zPt8jz4aTjopdjXpoX17uOUWmDzZl7HPUQpgkcYUgo942LjRVwrOhtUuUmXYMD8pN3IkLFkSu5ooFMAijenRR2HiRLjpJthjj9jVpJe8PHjgAf/lNGKE/7LKMQpgkcayeDFcfDEcfLDfypZ69fJREc8952GcYxTAIo1h0yYYOtRHPzz4IOTnx64ofV10ka8EfdFFvjBpDlEAizSG666DN9/0q97U9bBteXkwbhy0awc/+lFOLV+kABZJtfHjfa7fc8/1VrDUrFs3eOwx+PBD+MlPcmZomgJYJJUmT4b//E8YONDH/ErtHXmkzxUxfrz/BZEDCmIXIJI15szxVYB79oQJE6Bly9gVZZ7LL4ePPvK5InbeGc47L3ZFjUoBLJIKH37oLTgzP6PfqVPsijKTGYwe7SNIzj/fT16OGBG7qkajLgiRhpoxw8/il5fDa6/ppFtDFRTAU0/5/MHnnOMXsGQpBbBIQzz/vPf3msGrr8I++8SuKDu0aOEzxw0e7MPTRo6EsrLYVaWcAlikPsrK/Oq2E07wiwnefhv+7d9iV5VdKkL40kth1Cg47jhfzimLKIBF6mrBAu9yuO46H7f6+uvQo0fsqrJTfr6PJhkzxj/nPn3g2WdjV5UyCmCR2lq71mc122cfn8f24Yfh8cehTZvYlWW/4cNh2jT/Rfcf/+H9w3Pnxq6qwRTAIjVZtcqHRe26q89j+8MfwqxZPt5Xs5s1nX32galT4fbb4Y03/P6ZZ/q/RYZSAItUp7QUXnkFzjrLr9K65ho48EC/vPixx9TlEEvz5j5WeO5cuOQS7yPed1847DBf7n7lytgV1omFHJwCTrZUVFQUiouLY5cRT1mZX0jx5ps+muHFF+HLL7174cc/9uFQ/frFrlI2t2yZz6I2dizMm+d9xocdBoMG+W2/ftC6daO9vZlNCyEU1fv5CmCBHAjgjRu9dbRsGXzxhQ/0X7QI5s/34J0xA9at82O7doWjjvK+xmOOadQfYEmREOCdd+CZZ+Bvf/Pln8C7iPbc00eo7LEHFBb6lYrduvnCqB07+r9vPbuSFMCSEkWdO4fiwYMb/4229v+tYn/Vx0PYcisv99uyssqttNS3TZtgwwbf1q/3bc0a+Ppr/7o6PXrAXntB797eWurf339g1beb2Vas8Hk5pk/3E6Zz5vgv240btzy2oADatvWtVSvYbjsfAteihXd5NGvmx+TnV2477QR33KEAltQoat48FDfVUulbC7eK/VUfN/v2lpfnt1V/GPLz/YekWbPKH5zttvOtbVvvRmjf3rfOnaFLFw/eHj38WMkN5eX+188nn/htSYn/VbRypf+S/vpr/yto/frKX+QbN/ov9tLSyl/45eW+uvWLLzY4gDUXhLi+fSGbuyBE8vK866Fbt9iV/ItGQYiIRKIAFhGJRH3AAoCZfQ3MiV1HDToBy2MXUYNMqBEyo85MqHGvEELb+j5ZfcBSYU5DTiY0BTMrVo2pkQl1ZkqNDXm+uiBERCJRAIuIRKIAlgpjYhdQC6oxdTKhzqyvUSfhREQiUQtYRCQSBbCISCQK4BxnZsea2Rwzm29mV8eup4KZ9TSzV81slpnNNLNLkv0dzewlM5uX3HZIg1rzzWy6mU1M7u9qZlOTz/QJM2seub72ZvaUmX1oZrPN7KB0+xzNbGTy7zzDzB43s5bp8Dma2QNmVmJmM6rsq/azM/c/Sb3vm1mN85cqgHOYmeUDdwM/APYBTjOzdFnWtxS4PISwDzAAuCCp7WpgUghhD2BScj+2S4DZVe7fAtwZQugFrATOjlJVpbuAF0IIewP74bWmzedoZt2Bi4GiEEJvIB8YSnp8jn8Ejt1s39Y+ux8AeyTbCODeGl89hKAtRzfgIODFKvevAa6JXddWap0AfB+/Wq9bsq8bfgFJzLp6JD+ERwITAcOv3iqo7jOOUN/2wEKSE+5V9qfN5wh0Bz4FOuIXh00EjkmXzxEoBGbU9NkB9wGnVXfc1ja1gHNbxX/8CouTfWnFzAqBA4CpQJcQwufJQ0uBLpHKqjAKuBIoT+7vAHwVQihN7sf+THcFlgEPJt0k95tZa9LocwwhLAFuBz4BPgdWAdNIr8+xqq19dnX+eVIAS1ozszbA08ClIYTVVR8L3syINo7SzE4ASkII02LVUAsFQD/g3hDCAcBaNutuSIPPsQMwGP9lsRPQmi3/7E9LDf3sFMC5bQnQs8r9Hsm+tGBmzfDwHRdC+HOy+wsz65Y83g0oiVUfcDDw72b2MfAnvBviLqC9mVXMsxL7M10MLA4hTE3uP4UHcjp9jkcBC0MIy0IIm4A/459tOn2OVW3ts6vzz5MCOLe9A+yRnG1ujp/4+EvkmgA/owyMBWaHEO6o8tBfgGHJ18PwvuEoQgjXhBB6hBAK8c/ulRDC6cCrwMnJYbFrXAp8amZ7JbsGAbNIo88R73oYYGatkn/3ihrT5nPczNY+u78AZyajIQYAq6p0VVQvVse7tvTYgOOAucBHwHWx66lS1yH4n3bvA+8m23F4H+skYB7wMtAxdq1JvYcDE5OvdwPeBuYD44EWkWvbHyhOPstngQ7p9jkCvwI+BGYAjwAt0uFzBB7H+6U34X9NnL21zw4/AXt38rP0AT6qY5uv36BLkc3sWPxPrnzg/hDCzfV+MRGRHFPvAE7GkM7FhwYtxv+cPS2EMCt15YmIZK+GTMj+XWB+CGEBgJn9CT+TudUA7tSpUygsLGzAW0pDffWVL/7as+e390+bNm15CKFzxf3v5/1IszSJVPFS+fitLOddfw0J4OrGvPXf/CAzG4FfFcLOO+9MsVbejerKK2H06C0XQDazRXEqEsldjT4KIoQwJoRQFEIo6ty5c81PkEa1aRM0axa7ChGBhgVwWo8hlept3AjNo04NIyIVGhLAaTuGVLZu3Tpo3Tp2FSICDegDDiGUmtmFwIv4MLQHQggzU1aZNIrVq6FNm9hViAg0cFn6EMJzwHMpqkWawIoVsMMOsasQEdClyDlnyRLo2jV2FSICCuCcUlYGH38Mu+8eu5KIzHwTSQMK4BwyYwaUlkLv3rErERFoYB+wZJa//91vDzssbh1RVVx6X9EKtqQNEsq//bhIE1ALOIc8+ijss8+WlyGLSBxqAeeIKVPgnXf8MmShsqUbyvw2Lx8Aa+4/EmFTshJOeVlTVyY5RC3gHLBxI5x3HnTpAmeeGbsaEamgFnCWCwGuuw7efReefRbatYtdUZpKWrphg99agf9o2HZ+2WDYuMlvN22MUJxkK7WAs1gIPvvZ7bd7C3jw4NgViUhVagFnqVWr4OKL4eGH4YIL4H/+J3ZFmSWUln7rNq9lS7/d3i8jDOvWA1C+bl2E6iRbqAWchZ591kc7PPooXH+9n3jL07+0SNpRCziLvPUW/Pa38Ne/Qt++HsTf+U7sqrJD+Tff+BfJbX7SmV6wi4/pC6vXAFC2cmXTFycZS+2iDFdWBs88AwcfDAMHwj/+ATfd5CteKHxF0ptawBlq9mx44gkYNw7mz4fCQu/n/a//0nSTTaFs9Wr/IrnNT1Z7yeu7t98u+wqA0s+XNn1xkjEUwBlk7lwP3Sef9HkdzPyy4htvhB/+EAr0rymSUfQjm8bWrIE33oCXX4aXXoIPPvDQPeQQP7E2ZAh06xa7SgEoW7bMv0hu8wp39v1H9AOg+WJvEZfNW9D0xUnaUgCnkU2bYOpUD9xJk/zy4dJSX8Pt4INh1Cg4+WTo3j12pSKSCgrgiEpKPHArtrfegrVrvZV74IFw+eVw1FEevtttF7taqYvSjz8BID+5Zd+9APh66AAA2s37GoAwTat45TIFcBPZsMEvB54yxcN2yhRYuNAfy8+HPn1g2DAYNAgOPxw6doxarog0AQVwI1izBt5/3wN3+nS/ff99nxQHvAthwAC/PLh/f2/taqXi7FY2cw4AbZMGb/nA/QD4YuRAADq9twGAglemNX1xEo0CuIG++OLbQTt9OsybVznbYceOcMABcMklHrb9+0OPHnFrFpH0oACupfXrYdYsH4lQdVtaZZhnYaGH7emnw/77+9c9emgJMtmSTX4PgK6T/f6G4/yqmUW3HQRAt7d8hY5Wf57a9MVJk1EAb6asDBYs2DJo58+H8mTVmpYtfa6FY46pDNr99oP27ePWLiKZJacDuKRky6CdORMqJrgy8xWE+/SBoUP9tk8f6NXLT5yJpEqL594BYPfn/P7qH/toiYV/6gvA9v/nJwk6PvBW0xcnjSYnAnjdOg/WzcO2pKTymM6dPVyHD68M2n331ckxEWk8WRfAX37pJ8L++c/K27lzK0+KtWzpy7Iff3xl0Pbp48v1iKSLdo9NSW79fsn5Plqi9es+58RHz+wBQNc7Jzd9cZIyGRvAIcBnn20Ztp98UnlMz57Qr9+3uw92313dByKSHjImgMvLvRvh9dd9foTXX4fPP698fM894aCDfPWHAw7wrVOnePWKpNKO93hLd+09fr/sam8BnzDT5x++90/HA9DzN2oRZ5K0DeBNm2DatMqwffNNqJjrunt3v1pswABv4e63H7RtG7VcEZE6S6sADsEnFH/4YZ9ysWLK1T339OkWDz3Up18sLNTYWslt3W/2lu7EmzsAUPYbP8lxy0IfN3zKE5cCsOvVGjWRzmoMYDPrCTwMdAECMCaEcJeZdQSeAAqBj4FTQgj1Wo9lwQJ45BEP3gULfOTBkCFw4ok+9WLXrvV5VRGR9FabFnApcHkI4Z9m1haYZmYvAWcBk0IIN5vZ1cDVwFV1LeB3v4MrrvAW7aBBcMMN3trV8C+R2iv8hbd0r/pFfwDKb/P9Ty/20RT7Jy3i3S+f0vTFyVbVGMAhhM+Bz5Ovvzaz2UB3YDBweHLYQ8Br1DGAx4zx8B0yBO6800ctiIjkijr1AZtZIXAAMBXokoQzwFK8i6LWFi6Ec8+FoiJ47DGfdFxEUmP3n3mLeMjP/Iq6cKfvf/Gzd/3xP50LQK/L1CKOqdarIptZG+Bp4NIQwuqqj4UQAt4/XN3zRphZsZkVL6tYtgWfpOaII3z87nPP1a94EZFMVqsANrNmePiOCyH8Odn9hZl1Sx7vBpRU99wQwpgQQlEIoahzsnIsQLNm8Oyz3gIeMgSOO84XnPzmmwZ9PyJSjV4jp9Br5BSO2Wl/jtlpfyyABe8jfnrxFD667SA+SmZik6ZTYwCbmQFjgdkhhDuqPPQXYFjy9TBgQl3fvG1beOEFuOoqn5th6FAf8TBihI/7rZh9TEQkG1kI1fYcVB5gdgjwBvABUBGJ1+L9wE8COwOL8GFoX27rtYqKikJxcXG1j5WVwWuvwUMPwdNP+wQ6HTr4MLTDDvMxwP36ectZUs/MpoUQiirufz/vR9v+jyFZZeHN3vp98tRRAJz6qI+aqBhdIfBS+fiUX31Qm1EQ/wC29saDUlVIfr4PQxs0CO65ByZMgFdf9Svh/vpXP6ZVK7/cuOKCjO9+V8PVRCRz1dgCTqVttYC3ZelSv0KuYh6I997zq+bMYK+9fN6Hfv0q54DQgpZ1pxawVPXJ9T772vmn/g2A+x7xuSYqrsDLRVFawOmga1c4+WTfAL76CiZPhnfe8VEU//gHPP545fG77FIZyP36+aoVO+2ky5dFJL1kRABvrn17HzVx3HGV+5Yv9zCu2P75Tx9lUdHA79ChckrKvn39tndvTeIjUp2df53MNfFrn2sif6Tvr5iP+OPHewHQ+V71ETdERnRB1NfXX3t3xXvv+SiL99+HGTN8f4XCwm9PzN6nj0/+k2sn+9QFIXXx5U/8pN2qo9cCsMNfWgGVE8lno5ztgqivtm19FMUhh1TuCwEWLdpyeaLnn4fSUj+meXPYe+8tg1krHItIKmV1C7guNmyAOXO8lVw1mBcvrjymfXvvtti8K6Ndu3h1p4pawNIQa4f4JEBLB/ilBT1e9dZMxWKj2UAt4EbUooUHat++396/cqV3W1QN5XHjKucqBl8luWJ5+orbrl3VWhaRbVMLuB5CgE8/9dbyu+/6Nn26z2VcYccdKwN5//3hwAM9qNM1lNUCllQqPfJAAJbv14IuU72f2Ca/F7OkBlMLOE2Ywc47+3bCCZX7V63yUJ4+vTKU77jDl1cCH5/cv79vAwb4hSQdOsT5HkQkPgVwCm2/vV+ld+ihlfs2boRZs6C4GKZOhSlTfP6Lij889tzTw7h/fxg40LtA8mo9R51Ieip4ZRoAXV8BO3BfAFYP9akx23/gC+eUzZwTp7g0ogBuZM2bV3ZD/PSnvm/1ag/kKVM8lF94wZdjAl/J+Ygj4Kij/LLs3XZL324LEWkYBXAE7drBkUf6BpVD415/HSZN8m38eH+ssLByjoxjj1WXhWSeMG0mAG2nJTv22A2AsiP6AdBigc8TXrro0yavLTYFcBow86AtLIQzz/RAnjOnMoyffhrGjvWLQ44+Gk45BQYP9i4PEclcCuA0ZOYXguy9N1xwgU/V+c47HsRPPgl/+5t3bRx7rIfxkCHQsmXsqkVqp2yeDxfKn5fs6L4TAHl99wbAPl/hx1VZQSdb6XRPBsjP9xN1t90GH38Mb73lwTxtGpxxhk8+9JvfwIoVsSsVkbpQAGcYMw/jO+6ATz6Bl1/2McbXX++rSl94ofcni2SK0iWfUbrkM8rf/5Dy9z/0OQFKSynYpScFu/Qkr21b8rJ01iwFcAbLy/OTc889V7mk05gxsM8+MGqUd12ISPpSAGeJ3r3hgQdg3jw4/HAYOdLHFc+cGbsykbopW7mSspUrKV30qY+MKC+H8nLyO+1AfqcdyGvZkrwsOemhAM4yu+wCEyfCY4/5pdEDB/oCpyKSfhTAWcgMTjvNJ6Xv2tWHrk2aFLsqkfopX7uW8rVrKVu+grLlKwghEEIgr3Vr8lq3xgoKsILMHNClAM5iPXv6xR277ebLOS1ZErsiEalKAZz3yuHQAAAGyklEQVTlunSBZ57x+Y6HD6+cg0IkU4UNGwgbNvyrZVzBWrTAWrTwPwEz5Pp9BXAO6NULbrrJV/149dXY1YhIBQVwjjj3XJ8O8557YlciklqhtNS3pGX8L3n5vqUxBXCOaNkShg2DCROgyl9tIhKRAjiHHH20X2Q0JXsXrhXxEx0hQHmZbxV9wmnYN6wAziEH+UrivP123DpExGXm4Dmpl+2391ERH30UuxKRJpTGQ3/UAs4xhYU+iY+IxKcAzjGdOmnaSpF0oQDOMR06wMqVsasQEahDAJtZvplNN7OJyf1dzWyqmc03syfMrHnjlSmp0qoVrF8fuwoRgbq1gC8BZle5fwtwZwihF7ASODuVhUnjaNlSASySLmoVwGbWAzgeuD+5b8CRwFPJIQ8BJzVGgZJazZrBpk2xqxARqH0LeBRwJVCe3N8B+CqEUJrcXwx0r+6JZjbCzIrNrHhZDiyyl+4KCvxiDBGJr8YANrMTgJIQwrT6vEEIYUwIoSiEUNS5c+f6vISkUF6eLzAgIvHV5kKMg4F/N7PjgJZAO+AuoL2ZFSSt4B6AZpvNAApgkfRRYws4hHBNCKFHCKEQGAq8EkI4HXgVODk5bBgwodGqlJQxS+sLg0RySkPGAV8FXGZm8/E+4bGpKUkakwJYJH3UaS6IEMJrwGvJ1wuA76a+JGlMaTYZlEhO05VwIiKRKIBzjFrAIulDASwiEokCWEQkEgWwiEgkCmARkUgUwCIikSiARUQiUQCLiESiABYRiUQBLCISiQJYRCQSBbCISCQKYBGRSBTAIiKRKIBFRCJRAIuIRKIAFhGJRAEsIhKJAlhEJBIFsIhIJApgEZFIFMAiIpEogEVEIlEAi4hEogAWEYlEASwiEokCWEQkEgWwiEgkCmARkUgUwCIikdQqgM2svZk9ZWYfmtlsMzvIzDqa2UtmNi+57dDYxYqIZJPatoDvAl4IIewN7AfMBq4GJoUQ9gAmJfdFRKSWagxgM9seOAwYCxBC2BhC+AoYDDyUHPYQcFJjFSkiko1q0wLeFVgGPGhm083sfjNrDXQJIXyeHLMU6FLdk81shJkVm1nxsmXLUlO1iEgWqE0AFwD9gHtDCAcAa9msuyGEEIBQ3ZNDCGNCCEUhhKLOnTs3tF4RkaxRmwBeDCwOIUxN7j+FB/IXZtYNILktaZwSRUSyU40BHEJYCnxqZnsluwYBs4C/AMOSfcOACY1SoYhIliqo5XEXAePMrDmwAPgvPLyfNLOzgUXAKY1ToohIdqpVAIcQ3gWKqnloUGrLERHJHboSTkQkEgWwiEgkCmARkUgUwCIikSiARUQiUQCLiESiABYRiUQBLCISiQJYRCQSBbCISCQKYBGRSBTAIiKRKIBFRCJRAIuIRKIAFhGJRAEsIhKJAlhEJBIFsIhIJApgEZFIFMAiIpEogEVEIlEAi4hEogAWEYlEASwiEokCWEQkEgWwiEgkCmARkUgUwCIikSiARUQiUQCLiERSqwA2s5FmNtPMZpjZ42bW0sx2NbOpZjbfzJ4ws+aNXayISDapMYDNrDtwMVAUQugN5ANDgVuAO0MIvYCVwNmNWaiISLapbRdEAbCdmRUArYDPgSOBp5LHHwJOSn15IiLZq8YADiEsAW4HPsGDdxUwDfgqhFCaHLYY6F7d881shJkVm1nxsmXLUlO1iEgWqE0XRAdgMLArsBPQGji2tm8QQhgTQigKIRR17ty53oWKiGSb2nRBHAUsDCEsCyFsAv4MHAy0T7okAHoASxqpRhGRrFSbAP4EGGBmrczMgEHALOBV4OTkmGHAhMYpUUQkO9WmD3gqfrLtn8AHyXPGAFcBl5nZfGAHYGwj1ikiknUKaj4EQgi/BH652e4FwHdTXpGISI7QlXAiIpEogEVEIlEAi4hEogAWEYlEASwiEokCWEQkEgWwiEgkCmARkUgUwCIikSiARUQiUQCLiESiABYRiUQBLCISiQJYRCQSBbCISCQKYBGRSBTAIiKRKIBFRCJRAIuIRKIAFhGJRAEsIhKJAlhEJBIFsIhIJApgEZFIFMAiIpEogEVEIlEAi4hEogAWEYlEASwiEokCWEQkEgWwiEgkCmARkUgUwDlmwAC45JLYVYgIgIUQmu7NzJYBi5rsDaUudgkhdI5dhEguadIAFhGRSuqCEBGJRAEsIhKJAlhEJBIFsIhIJApgEZFIFMAiIpEogEVEIlEAi4hEogAWEYlEASwiEokCWEQkEgWwiEgkCmARkUgUwCIikSiARUQiUQCLiESiABYRiUQBLCISiQJYRCQSBbCISCQKYBGRSBTAIiKR/D+iA6AGU5wC6AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x360 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#%% EMD\n",
    "\n",
    "G0 = ot.emd(a, b, M)\n",
    "\n",
    "pl.figure(3, figsize=(5, 5))\n",
    "ot.plot.plot1D_mat(a, b, G0, 'OT matrix G0')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solve Sinkhorn\n",
    "--------------\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "It.  |Err         \n",
      "-------------------\n",
      "    0|7.958844e-02|\n",
      "   10|5.921715e-03|\n",
      "   20|1.238266e-04|\n",
      "   30|2.469780e-06|\n",
      "   40|4.919966e-08|\n",
      "   50|9.800197e-10|\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFgCAYAAACFYaNMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XmcVNWZ//HP0900TbMjyK7t7ii4IAnuo+KuGTLRGI35iYlxj1tMNGom0cnERGNcxqiJ0RiN6KgYlxCXqNG4QgQxiiKIKLKDitAs0tv5/fHc29003XTT26nl+3696lVU1a1bDwX9radPnXuuhRAQEZHOVxC7ABGRfKUAFhGJRAEsIhKJAlhEJBIFsIhIJApgEZFIFMAikZnZAWY2qwP2e7KZ/a2F255qZi9v7mPSNgpgyXpJQLxtZmvNbImZ3WZmfZLHfmtmq5NLhZlV1rv9ZCfUFsxs+01tE0J4KYSwUyv3v7+ZvWpmK83sMzN7xcy+lOx3Qgjh8NbsVzqHAliympldDFwD/BDoDewNbA08Y2bFIYSzQgg9Qgg9gKuBB9LbIYSj4lXuzKyoDc/tBUwCbgb6AUOBq4D17VNd+2vL3zcXKYAlayUBdBVwXgjhqRBCZQjhI+AEoAz4Viv2eZCZLTCzS8xsmZktNrOvmtnRZjY76TIvr7f9l83sNTP7PNn2N2ZWnDz2YrLZv5KO+xv19n+pmS0B7krvS56zXfIao5LbQ8xsuZkd1Ei5OwKEEO4PIVSHENaFEP4WQngree4GQwdJN36Wmb2f1HuLmVkT78OvzOxlM+td777rzGyFmX1oZkfVu3+ImT2e1D3HzE6v99iVZjbRzO41s1XAqcl9D5rZPWZWbmbvmNnozfynygkKYMlm+wIlwJ/r3xlCWA08ARzWyv0OSvY7FPgJ8Hs8zPcCDgD+y8y2SbatBi4C+gP7AGOBc5I6Dky22T3puB+ot/9+eKd+RoPaPwAuBe41s1LgLuDuEMILjdQ5G6g2s7vN7Cgz69uCv9uxwJeA3fAPqiPqP2hmBWb2++Txw0MIK5OHxgCzkr/ntcCd9cL7/4AFwBDgeOBqMzuk3m7HAROBPsCE5L7/SJ7XB3gc+E0Las85CmDJZv2BT0IIVY08tjh5vDUqgZ+HECrxkOgP3BRCKA8hvAO8C+wOEEKYFkKYHEKoSrrv3wH/3sz+a4CfhhDWhxDWNXwwhPB7YA4wBRgMXNHYTkIIq4D9gYB/SCxPOtGBm3jtX4YQPg8hfAw8D+xR77EuwP34h8NXQghr6z02L4Tw+xBCNXB3UtdAMxsO7AdcGkL4IoTwJnAHcEq9574WQng0hFBT7+/7cgjhiWR/fyJ5P/ONAliy2SdA/ybGFQcnj7fGp0kwAKSBsbTe4+uAHgBmtqOZTUq+/FuFjzM3F/zLQwhfNLPN74ERwM0hhCbHdEMIM0MIp4YQhiXbDwFu3MR+l9T789r075HYHu9WrwohVDT1vHrB3CN5vc9CCOX1tp2H//aQmt+COkrycXxYASzZ7DX8C6ev1b/TzHoARwHPdUINtwHvATuEEHoBlwONjqvWs8klCJP6bwTuBK40s34tKSSE8B7wRzyIW2Mm8G3gSTNr6ayMRUA/M+tZ776tgIX1S2tlPTlPASxZKxmfvAq42cyONLMuZlYGPIiPSf6pE8roCawCVpvZzsDZDR5fCmy7mfu8CZgaQvgu8Ffgt41tZGY7m9nFZjYsuT0cOAmYvJmvVyuEcD/+IfKsmW3Xgu3nA68CvzCzEjPbDTgNuLe1NeQTBbBktRDCtXhgXIcH4RT8V96xm/rVvR39APgmUI4PGzzQ4PErgbuTWQcnNLczMxsHHEldkH8fGGVmJzeyeTn+5dgUM1uDB+8M4OJW/D1qhRDuBv4b+Hvygdack/BZJ4uAR/Dx7WfbUkO+MC3ILiIShzpgEZFIFMAiIpEogEVEIlEAi4hEkncTn6Vx/fv3D2VlZbHLEMkq06ZN+ySEMKC1z1cACwBlZWVMnTo1dhkiWcXM5rXl+RqCEBGJRAEskk9WroRXX4V580DHAESnABbJdSHAQw/BiBHQpw/stx+UlcGWW8KPfgTl5c3uQjqGAlgkly1ZAmPHwgknQGEh/Pzn8NhjcOutcPDBcM01sNNO8PTTsSvNS/oSTiRXLVjg4btgAdxyC5x5podw6uyzYcoUOP10+MpX4P774bjj4tWbh9QBi+SihQvhgAO8A/7b3+CcczYM39SYMfDSS/ClL3mX/NBDnV9rHlMAi+SaykoP008+geee8zHfTend24cg9tkHxo+HGTM6p05RAIvknEsu8ZkOd94Jo1t4rssePWDiRA/j446DVas6tkYBFMAiuWXSJLjxRjj/fO+CN8egQfDAA/DBB3DuuR1Tn2xAASySK1av9rHeESPgV79q3T4OPBAuvxzuvReeeaZ965ONKIBFcsWVV8L8+fC730Fxcev3c/nlsMMOHuZfNHfuUGkLBbBILnjrLR96OP102Hfftu2rpMTnCc+ZA7/4RfvUJ41SAIvkgksugV694Je/bJ/9HXoofOMbPpSxaFH77FM2ogAWyXbPP+/TyC6/HPq16Az2LXP11VBVBf/93+23T9mAAlgkm4Xg6zkMGwbf+1777nvbbf3ouTvugNmz23ffAiiARbLbI4/AP/8JV13lY7ft7cc/9v3+13+1/75FASyStUKA//kfn7Fwyikd8xoDB8KFF/ohyjNndsxr5DEFsEi2evJJmD4dLrsMijpwXa0LL4Ru3drvCz6ppQAWyUZp97vVVvCtb3Xsa/XvD2ecARMmwNy5HftaeUYBLJKN/vEPeO01n37WpUvHv94PfuCrqV17bce/Vh5RAItko2uv9TNafOc7nfN6Q4f6Sml//CMsW9Y5r5kHFMAi2ebdd33893vf87HZzvL978P69X6UnLQLBbBItrnxRp8adtZZnfu6O+8MxxzjAbxuXee+do5SAItkk2XL4J57fNrZgAGd//oXXwzLl/sXctJmCmCRbPLb3/owwEUXxXn9gw6CPfaA66/Xae3bgQJYJFtUVMBtt8FRR/lwQAxmHv4zZ/rpjqRNFMAi2WLiRD/J5nnnxa3jG9/w4Y+bb45bRw5QAItki5tv9sOOjzgibh1du/qBGX/5C3z4YdxaspwCWCQbTJ0Kkyf71LOCDPixPftsr0NT0tokA/4lRaRZN9/sZy4+9dTYlbihQ/3syXfcAWvXxq4maymARTLdJ5/42YpPOcXPepEpvvc9+PxzuP/+2JVkLQWwSKa7806fenbOObEr2dD++8PIkXDLLZqS1koKYJFMVl3tU88OOgh23TV2NRsy8w+F6dN9fFo2mwJYJJM98QTMmwfnnhu7ksZ961s+LPKb38SuJCspgEUy2S23wJAhMG5c7Eoa16OHr5L20EOwdGnsarKOAlgkU73/vp/t+MwzO2fN39Y65xyorPQZEbJZFMAimeq22/xUQ6efHruSTdt5Zxg7Fn73Oz+NvbSYAlgkE61dC3fd5XNtBw+OXU3zzj0X5s+HSZNiV5JVFMAimei++3yObaZ++dbQV74Cw4f7mLW0mAJYJNOE4EE2cqTPtc0GRUW+QPyzz8J778WuJmsogEUyzSuvwJtv+pFmZrGrabnvfheKizUlbTMogEUyzc03Q58+cPLJsSvZPFtuCSeeCHffDatWxa4mKyiARTLJwoXw8MNw2mnQvXvsajbfeefB6tV+9mRplgJYJJP89rdQU5N56z601OjRMGaMD0PU1MSuJuMpgEUyxRdf+FzaY4+FbbeNXU3rnXeeH0Ty1FOxK8l4CmCRTDFhgp9x+MILY1fSNl//us9dvuGG2JVkPAWwSCYIwQNrt93g4INjV9M2xcU+g+PZZ+Htt2NXk9EUwCKZ4Nln4Z13/IzD2TT1rClnngndusFNN8WuJKMpgEUywQ03wMCBcNJJsStpH1ts4WfwuPdeWLYsdjUZSwEsEtvbb8OTT/phx127xq6m/Vx0EVRU6PT1m6AAFont2mt9zm+2rPvQUjvt5OsY33KLzw2WjSiARWKaN89Pann66dCvX+xq2t+ll8KKFfD738euJCMpgEViuv56/9Lt+9+PXUnH2Htv+Pd/979nRUXsajKOAlgklmXLvDM8+WRfyjFXXXopLFgAf/pT7EoyjgJYJJZf/cpPN3/ZZbEr6VhHHgl77QU//7mfukhqKYBFYli2DG69Fb75Tf+yKpeZwZVXwocfqgtuQAEsEsN11/naDz/+cexKOscxx6gLboQCWKSzLVniU7NOOin3u99U2gXPnaulKutRAIt0tquu8hkBV14Zu5LOdcwxPiviyiv9pKOiABbpVLNn+8yHM8+E7bePXU3nMvODThYt0hoRCQWwSGe6/HJfpOYnP4ldSRwHHOBnUP7lL+HTT2NXE50CWKSzvPiin27oBz/w86flq1/8wg9N/ulPY1cSnQJYpDNUVfkauVtvDT/8Yexq4tp1Vz/l0m23+dmf85gCWKQz3HKLr3p2ww1QWhq7mvh+9jNfsvLcc/P63HEKYJGOtmiRj/kefjh89auxq8kMffrANdfAq6/6aezzlAJYpCOF4DMeKir8TMG5cLaL9jJ+vH8pd9FFsHBh7GqiUACLdKR774VJk+Dqq2GHHWJXk1kKCuAPf/APpzPO8A+rPKMAFukoCxbA+efDfvv5tWxs++19VsQTT3gY5xkFsEhHqKyEE0/02Q933QWFhbErylznnedngj7vPD8xaR5RAIt0hCuugFde8aPeNPSwaQUFMGEC9OoFX/96Xp2+SAEs0t4eesjX+j3rLO+CpXmDB8N998F778F3vpM3U9MUwCLt6dVX4f/9P9h3X5/zKy13yCG+VsRDD/lvEHmgKHYBIjlj1iw/C/Dw4fDYY1BSErui7HPxxfDBB75WxFZbwdlnx66oQymARdrDe+95B2fm3+j37x+7ouxkBjff7DNIzjnHv7w844zYVXUYDUGItNWMGf4tfk0NvPCCvnRrq6IimDjR1w8+80w/gCVHKYBF2uLJJ3281wyefx522SV2Rbmha1dfOW7cOJ+edtFFUF0du6p2pwAWaY3qaj+67dhj/WCCf/4T/u3fYleVW9IQvvBCuPFGOPpoP51TDlEAi2yuuXN9yOGKK3ze6osvwrBhsavKTYWFPpvk9tv9fR45Eh59NHZV7UYBLNJSa9b4qma77OLr2N5zD9x/P/ToEbuy3Hf66TBtmn/Q/ed/+vjw7Nmxq2ozBbBIc1au9GlR22zj69h+7Wvw7rs+31erm3WeXXaBKVPguuvgpZf89imn+L9FllIAizSmqgr+/nc49VQ/Suuyy2Cvvfzw4vvu05BDLMXFPld49my44AIfI951VzjwQD/d/YoVsSvcLBbycAk42djo0aPD1KlTY5cRT3W1H0jxyis+m+Hpp+Gzz3x44Zvf9OlQo0bFrlIaWr7cV1G78054/30fMz7wQBg71q9HjYLu3Tvs5c1sWghhdKufrwAWyIMArqjw7mj5cli61Cf6z5sHc+Z48M6YAWvX+raDBsGhh/pY4xFHdOgPsLSTEOD11+GRR+Cvf/XTP4EPEe24o89Q2WEHKCvzIxUHD/YTo/br5/++rRxKUgBLuxg9YECYOm5cx79QU//f0vvrPx7CxpeaGr+urq67VFX5pbIS1q/3y7p1flm9GsrL/c+NGTYMdtoJRozwbmnMGP+B1dhudvv0U1+XY/p0/8J01iz/sK2o2HjboiLo2dMvpaXQrZtPgeva1Yc8unTxbQoL6y5DhsD11yuApX2MLi4OUzvrVOlNhVt6f/3HzTa8FBT4df0fhsJC/yHp0qXuB6dbN7/07OnDCH36+GXAABg40IN32DDfVvJDTY3/9vPxx369bJn/VrRihX9Il5f7b0Hr1tV9kFdU+Ad7VVXdB35NjZ/d+umn2xzAWgtC3G67QS4PQYgUFPjQw+DBsSuppVkQIiKRKIBFRCLRGLAAYGblwKzYdTSjP/BJ7CKakQ01QnbUmQ017hRC6NnaJ2sMWFKz2vJlQmcws6mqsX1kQ53ZUmNbnq8hCBGRSBTAIiKRKIAldXvsAlpANbafbKgz52vUl3AiIpGoAxYRiUQBLCISiQI4z5nZkWY2y8zmmNmPYteTMrPhZva8mb1rZu+Y2QXJ/f3M7Bkzez+57psBtRaa2XQzm5Tc3sbMpiTv6QNmVhy5vj5mNtHM3jOzmWa2T6a9j2Z2UfLvPMPM7jezkkx4H83sD2a2zMxm1Luv0ffO3P8m9b5lZs2uX6oAzmNmVgjcAhwF7AKcZGaZclrfKuDiEMIuwN7AuUltPwKeCyHsADyX3I7tAmBmvdvXADeEELYHVgCnRamqzk3AUyGEnYHd8Voz5n00s6HA+cDoEMIIoBA4kcx4H/8IHNngvqbeu6OAHZLLGcBtze49hKBLnl6AfYCn692+DLgsdl1N1PoYcBh+tN7g5L7B+AEkMesalvwQHgJMAgw/equosfc4Qn29gQ9JvnCvd3/GvI/AUGA+0A8/OGwScESmvI9AGTCjufcO+B1wUmPbNXVRB5zf0v/4qQXJfRnFzMqAPYEpwMAQwuLkoSXAwEhlpW4ELgFqkttbAJ+HEKqS27Hf022A5cBdyTDJHWbWnQx6H0MIC4HrgI+BxcBKYBqZ9T7W19R7t9k/TwpgyWhm1gN4GLgwhLCq/mPB24xo8yjN7FhgWQhhWqwaWqAIGAXcFkLYE1hDg+GGDHgf+wLj8A+LIUB3Nv61PyO19b1TAOe3hcDwereHJfdlBDPrgofvhBDCn5O7l5rZ4OTxwcCyWPUB+wH/YWYfAf+HD0PcBPQxs3Sdldjv6QJgQQhhSnJ7Ih7ImfQ+Hgp8GEJYHkKoBP6Mv7eZ9D7W19R7t9k/Twrg/PY6sEPybXMx/sXH45FrAvwbZeBOYGYI4fp6Dz0OjE/+PB4fG44ihHBZCGFYCKEMf+/+HkI4GXgeOD7ZLHaNS4D5ZrZTctdY4F0y6H3Ehx72NrPS5N89rTFj3scGmnrvHgdOSWZD7A2srDdU0bhYA++6ZMYFOBqYDXwAXBG7nnp17Y//avcW8GZyORofY30OeB94FugXu9ak3oOAScmftwX+CcwBHgK6Rq5tD2Bq8l4+CvTNtPcRuAp4D5gB/AnomgnvI3A/Pi5dif82cVpT7x3+Bewtyc/S2/isjk3uv02HIpvZkfivXIXAHSGEX7Z6ZyIieabVAZzMIZ2NTw1agP86e1II4d32K09EJHe1ZUH2LwNzQghzAczs//BvMpsM4P79+4eysrI2vKS01eef+8lfhw/f8P5p06Z9EkIYkN4+rODrWqVJpJ5nah5q4nTerdeWAG5sztuYhhuZ2Rn4USFstdVWTNWZd6O65BK4+eaNT4BsZvPiVCSSvzp8FkQI4fYQwugQwugBAwY0/wTpUJWV0KVL7CpEBNoWwBk9h1QaV1EBxVGXhhGRVFsCOGPnkErT1q6F7t1jVyEi0IYx4BBClZl9D3gan4b2hxDCO+1WmXSIVaugR4/YVYgItPG09CGEJ4An2qkW6QSffgpbbBG7ChEBHYqcdxYuhEGDYlchIqAAzivV1fDRR7DddrEricisdReRDqAAziMzZkBVFYwYEbsSEYE2jgFLdvnHP/z6wAPj1hFFU12stbAHafj0UNPgtg4clM2nDjiP3Hsv7LLLxochi0gc6oDzxOTJ8PrrfhiysFHnawXW6P0UNGh9a9JOt9Cvkk441DTogNUhSwuoA84DFRVw9tkwcCCcckrsakQkpQ44x4UAV1wBb74Jjz4KvXrFriiStANNx4LTDrWJjtfS7QoKNnxeEyzdf03SEVfXbPA6tR2yOmOpRx1wDgvBVz+77jrvgMeNi12RiNSnDjhHrVwJ558P99wD554L//u/sSvKEA074YbSTrXQNtjOCpMx3/S6YMPHG+7fqqv9ZtIJ1932a9LbDTtjdcR5RR1wDnr0UZ/tcO+98JOf+BdvBfqXFsk46oBzyGuvwS9+AX/5C+y2mwfxl74Uu6oMVdtppmO0/gllBcntkHS+6fZpx1uYbFeU/Oik14WNf8JZOhZcVeX7Ta5rb1c2uF+dcV5RX5TlqqvhkUdgv/1g333h5Zfh6qv9jBcKX5HMpg44S82cCQ88ABMmwJw5UFbm47zf/raWm9wsTXXCJGO1yRhvSMZwrLBBJ5p0vqE4Oc1Il6INtq9tcdJJF+kYcNL5Fqyv9O0rKvz+9ev9dlOdcfr8jeqXbKQAziKzZ3voPvigr+tg5ocV//zn8LWv1f02LCLZQT+yGWz1anjpJXj2WXjmGXj7bQ/d/ff3L9aOOw4GD45dZY5oohNOx2rTseDafjOdHdHgUy8U+fNqSrwjrin2WROhS4Mj75IxXqvw1yv4IumI165PbiedcHJd2xlXJB1z2iGnHbHGirOSAjiDVFbClCkeuM8954cPV1X5Odz22w9uvBGOPx6GDo1dqYi0BwVwRMuWeeCml9degzVrvLnaay+4+GI49FAP327dYlebJ9IOMiRjrqHxeb4N1R45V+Qdr3X1DjjtfKtK/f6qkqRDLtpwvwVVvt+iL0r9eq2/ftFq73gLVn+RXK/z/a5Lrtf72HHtGHLDsWJ1xBlNAdxJ1q/3w4EnT/awnTwZPvzQHysshJEjYfx4GDsWDjoI+vWLWq6IdAIFcAdYvRreessDd/p0v37rLV8UB3wIYe+9/fDgMWO829WZijNU0kHWdpTp/NyaDdd+qD0CLnlaQdIRFyRjwpR4B5x2vhU9/Lq6myX3p6/n2xVW+B1Fa7sC0LXcO+Piz33st8vKpCNe5Z1wweq1/vS0M07HkNP5xppXnJEUwG20dOmGQTt9Orz/ft3/7379YM894YILPGzHjIFhw+LWLCKZQQHcQuvWwbvv+kyE+pclS+q2KSvzsD35ZNhjD//zsGE6pVhOaDg2nN6/vkFHmc5uSDrjooadZkExANVd/UevMpmzXdnTr6u6+/Y1yZITljTeRV/4HV3K/brrCt9Pt8/8V6eun/pYcZfPvBMuWLXGy1qbdsbeMWv2RGZRADdQXQ1z524ctHPm1P22WVLiay0ccURd0O6+O/TpE7d2EckueR3Ay5ZtHLTvvANJ04CZn0F45Eg48US/HjkStt++blEsyVM1G86SSMdc0yPdLJ0/nFwXpUe+VSUtb02JP8/8R7Cm2PdTmXwXUN3Lty/s4dfVhb7fiir/j1e+2mdZdFnht0uW+35Kl/mYcbfl3lJ3/cTHhAs/W+2vl44Vp51xOq9YHXEUeRHAa9d6sDYM22XL6rYZMMDD9fTT64J211315ZiIdJycC+DPPvMvwt54o+569uy6D/SSEj8t+zHH1AXtyJF+uh6RzdZwbDg5so0G6/5a0mkWJms/dFvfI7ntn/CFld7RWnIE3tpkPnHo7h1wv17esW7Z3TvZomRweGWFTxBfuso73qXLfH8li32MuHSx77fH4uT+5T4WXPip78dWJZ3x2mT2RDJVJx0rVkfcsbI2gEOARYs2DtuPP67bZvhwGDVqw+GD7bbT8IGIZIasCeCaGh9GePFFXx/hxRdh8eK6x3fcEfbZx8/+sOeefunfP169kqcazBuuHRuu3HBM2JJOs+s6vy5a6x1xl7U+Nly4zn8011T67eVJA9q92Lcf2XcRAFv1+wCAgsHeqS4s6wvA258PAeCDJQMAWDXfO+XuC/x1eizy26WLvXMu+sQ74YKkI65Zk44RqyPuSBkbwJWVMG1aXdi+8gqsWOGPDR3qR4vtvbd3uLvvDj17Ri1XRGSzZVQAh+ALit9zjy+5uGqV37/jjr7c4gEH+PKLZWWaWysZbqM1JZLOcV1yXeljwZascla0zq97lCed8Cq/7rrKx3DLV3nHOnfNIADWbp3Mghjo+9+v5/sAHFA6G4Dj+vg42+wh/uXGS9vtBMCUxVsBsPij3v5683xsuMcC77S7L2rQEa8sB9QRd5RmA9jMhgP3AAPx+ee3hxBuMrN+wANAGfARcEIIYUVripg7F/70Jw/euXN95sFxx8FXvuJLLw4a1Jq9iohktpZ0wFXAxSGEN8ysJzDNzJ4BTgWeCyH80sx+BPwIuHRzC/j1r+EHP/COduxYuPJK73Y1/UtySjo2XLXhkWh1Y8JJR5ys5VCyxq+7rPSOtOQzXwui/FP/kf3s0y0BeHQbf3z+1j72e/QWbwGwX7ePABjRYwEAR3afB8CbW/jRQk8PHwnA89vuAMDiD/35PT/y1+k53+cTly7yH8S0I7akI66dNZGuOaF5xK3SbACHEBYDi5M/l5vZTGAoMA44KNnsbuAFNjOAb7/dw/e44+CGG3zWgohIvtisMWAzKwP2BKYAA5NwBliCD1G02IcfwllnwejRcN99vui4SN5obrZEcrsg6YRLVyZjwp95x1u63MeEy5d6xzp18Y4AvLud/xhOHTYHgGP7vAnAl7quBODw0srk9ssAHNH7bQAmDdkDgBe23h6AJR8kY8QfeQfcMxkjLk3mExd+4p2wrUznEW/iyDp1w01q8VmRzawH8DBwYQhhVf3HQgiBeuuTNHjeGWY21cymLl++vPb+YcPg4IN9/u4TT7SueBGRbGahBZ9OZtYFmAQ8HUK4PrlvFnBQCGGxmQ0GXggh7LSp/YwePTpMnTq19nZ5ORx2GLz+ui9sM348jBvnR6tJ5zKzaSGE0entwwq+rrYlhnR6j6VnYU7OsFHiY7JWmpwapbd3wlUD/HrNEH+8fLhvv3prH4vts41/L37wUJ8lkXa8e3T9HIBS8+2XVnvn/cq6MgCe+HQ3AKbN81kTBR/66/bwoWR6LvDtS5Z651v4aTI2XJ50xOnqaxUVOXN2jmdqHmr3uVfNdsDm51q5E5iZhm/icWB88ufxwGOb++I9e8JTT8Gll/raDCee6DMezjjD5/2mq4+JiOSiZjtgM9sfeAl4m/SUsXA5Pg78ILAVMA+fhvbZpvbVsAOur7oaXngB7r4bHn7YF9Dp29enoR14oM8BHjUKunTZnL+etJQ64AzVsCPu4l/bWPKlSUF3HwMOSUdcOcDHitcO3rAjXrOV/+h2L/Ox4DGD/Zj9A3rPAmDH4qUAdEnWmJiHDyNoAAAN7UlEQVRf5efEeqXcZ0m8vHRbAJbM2wKA0nleR8/5yX4X+dhv16XJbInPk454zZqNz86RpR1xR3TALZkF8TJ1Z1ppaGx7FVJY6NPQxo6FW2+Fxx6D55/3I+H+8hffprTUDzdOD8j48pc1XU1EsleLxoDby6Y64E1ZssSPkEvXgfjXv/zD0wx22snXfRg1qm4NCJ3QcvOpA84SLe2Ie3knnI4Rrx3kHfHqoUlHPDT55x3usyy2G+RfkI/s42tMDO3qY8c1wV/n4/X+Q/X2Cl9j4qPF3gkXpkfQzfe6ei7yLrfbknqrrjVccS1L5w5H6YAzwaBBcPzxfgH4/HN49VX/8m76dA/n+++v237rresCedQoP2vFkCE6fFlEMktWBHBDffrA0Uf7JfXJJx7G6eWNN+DRR+s+VPv2rVuScrfd/HrECC3iI1mmifWHa8/1lqzVYMnaDUXJkWu9lnlH3H2hX38xwDvX1YN9dsO8wVsDMGvwUN9uSz+n3MBe/vzexd69divysd4BW/j9y5OmpryLd9jV3ZJz3ZV6J15a2oXiT5LufEV6Vo7kfHXJOhi160tk6dhwW2RlADemf3+f0nbYYXX3lZf7cMW//uWzLN56y9ebKC+v26asbMOF2UeO9MV/9GWfZIVmAjk9oMOS0CtMvhzrkRzA0W2RB/L6/h7E6wZ4JKzd0g/EmD+gFwBz+/r+C3okC8sXJYdSJwFc3dNvr6tKTtFU4EMX1V2KKe3qfy7p6j9UBSuSxefLvSaSoQnycKGfnAngxvTs6bMo9t+/7r4QYN68jU9P9OSTkHxJS3Ex7LzzxsGsMxyLSHvK6QBujJl3vWVlvtpaav16mDXLu+Q0lP/xD5gwoW6bPn182KLhUEavXp39txBpQhPLYIYGp0aytenp670jLl2WdMS9kkOP+/nt9Vt4t7qun395t75vMsSQDN2FEn+99CwzaX9SnRxMtb6PYSGJmQK/s7jIO+Ki5EkFRf54WJMMTVhuTFtribwL4KZ07eqButtuG96/YgXMmLFhtzxhQt1axeBnSU5PT59eDxqkbllENk0B3Iy+fX3e8QEH1N0XAsyf793ym2/65Y03YOLEum223LIukPfYA/bay4NaoSydqrmOOPkiLP3SzpIvSIo/8zHhLkuSL9N6+e3K3t7FVvTx6Fjf07vZKt+M6uJ0mpxfFVRDdfJ9SkWP5LDq6g1X3ioqTKbUJePGte10Mn5NYwv81P+7ZTEFcCuYwVZb+eXYY+vuX7nSQ3n6dA/l6dPh+utrT35Av34wZoxf9t7bDyTp2zfO30FE4lMAt6PevTfulisq4N13YepUmDIFJk/29S/SD+8dd/QwHjMG9t3Xh0AKWrxGnchmatgRp81kgwXiSRbTSWdPFKzwrrWkm3fCXXskJ/ns4R1xVQ9vc6tKvXutKvH/xDVFdb/yFVQlS3AmDW5Nif+hpsqnsBUktdU+I/11Mf2BSGZJkC7hmQNjwwrgDlZcXDcM8d3v+n2rVnkgT57sofzUUz49Dnw63cEHw6GH+mHZ226rYQuRXKUAjqBXLzjkEL9A3dS4F1+E557zy0MP+WNlZXVrZBx5pIYspIM0OGVSbZfZcD5xMmfXypM5vV29e+2aLJdZXOodcU1JcXJdRCj2DjYUbNhJWGUyHp3cHbokrXHyXGuqs03H9NLSa9KuPvvGhhXAGaD+1LhTTvH/P7Nm1YXxww/DnXf6wSGHHw4nnODrJvfuHbtyEWkLBXAGMvMDQXbeGc4915uR11/3IH7wQfjrX31o48gjPYyPO06L2Es7a2r2RFUyjzgdj03HiosaLApU7B1yYXExITkCjmThoJDOcki/60gaV2u4AHg6OyLZd22Xnjxc20+n84VrGuwwCzphfd2TBQoL/Yu6X/0KPvoIXnvNg3naNPjWt3zxoZ/9DD79NHalIrI5smI5SmlcTY2vmfzrX/uh1N26wXe+Az/8oYfy5tBylNIqDZfHTMd501MpFRXVLplJ2iUXbXg77XQ3Up10smmHm45Pp7crG9yfjlt30JhwlFMSSeYqKPAv5554ou6UTrffDrvsAjfeWPv/UUQylAI4R4wYAX/4A7z/Phx0EFx0kc8rfued2JVJTgvBLzXVUFNNqKryS0UFoaKCmnVfULN6jV/KV/tlVblfVq6iZuUqQnpZvcYva9b65Ysv/FJZSais9I64ugZqgl8SZoaZeRduBViBeSee3MYsY+dyKoBzzNZbw6RJcN99MHeuh/Arr8SuSkQao1kQOcgMTjrJl+E89FCfuvb44z5cIdIp6s2iCOlE33SMtsF4ce3tdCw4PfKtua614dhuOv7c1NBbur8Mmh2hDjiHDR/uB3dsu62fzmnhwtgViUh9CuAcN3AgPPKIr3d8+ukZ9eEv+SIdJ244XlydXKoqCVWV1FQkly/WU/PFesL65FJRueElHWdOn19dQ2hkbLhWOhacgTKzKmlX228PV1/tU9Wefz52NSKSUgDnibPO8uUwb701diUiiSY646Y65NrblVUbXqqrN7zUhLq5wODzgUO9o+wyaFaEAjhPlJTA+PHw2GOQnPlFRCJTAOeRww/3g4gmT45diUgLNNMh13a2zV2a2m8GUADnkX328et//jNuHSLiNA84j/Tu7bMiPvggdiUi7SBDuti2UAecZ8rK4OOPY1chIqAAzjv9+2vZSpFMoQDOM337wooVsasQEdiMADazQjObbmaTktvbmNkUM5tjZg+YWXHHlSntpbQU1q2LXYWIwOZ1wBcAM+vdvga4IYSwPbACOK09C5OOUVKiABbJFC0KYDMbBhwD3JHcNuAQYGKyyd3AVzuiQGlfXbpsdFJZEYmkpR3wjcAl1J7tji2Az0MIyblAWAAMbeyJZnaGmU01s6nLly9vU7HSdkVFtWd0EZHImg1gMzsWWBZCmNaaFwgh3B5CGB1CGD1gwIDW7ELaUUGBn0tOROJryYEY+wH/YWZHAyVAL+AmoI+ZFSVd8DBAq81mAQWwSOZotgMOIVwWQhgWQigDTgT+HkI4GXgeOD7ZbDzwWIdVKe3GLCcOIBLJCW2ZB3wp8H0zm4OPCd/ZPiVJR1IAi2SOzVoLIoTwAvBC8ue5wJfbvyTpSBmyDKqIoCPhRESiUQDnGXXAIplDASwiEokCWEQkEgWwiEgkCmARkUgUwCIikSiARUQiUQCLiESiABYRiUQBLCISiQJYRCQSBbCISCQKYBGRSBTAIiKRKIBFRCJRAIuIRKIAFhGJRAEsIhKJAlhEJBIFsIhIJApgEZFIFMAiIpEogEVEIlEAi4hEogAWEYlEASwiEokCWEQkEgWwiEgkCmARkUgUwCIikbQogM2sj5lNNLP3zGymme1jZv3M7Bkzez+57tvRxYqI5JKWdsA3AU+FEHYGdgdmAj8Cngsh7AA8l9wWEZEWajaAzaw3cCBwJ0AIoSKE8DkwDrg72exu4KsdVaSISC5qSQe8DbAcuMvMppvZHWbWHRgYQlicbLMEGNjYk83sDDObamZTly9f3j5Vi4jkgJYEcBEwCrgthLAnsIYGww0hhACExp4cQrg9hDA6hDB6wIABba1XRCRntCSAFwALQghTktsT8UBeamaDAZLrZR1ToohIbmo2gEMIS4D5ZrZTctdY4F3gcWB8ct944LEOqVBEJEcVtXC784AJZlYMzAW+jYf3g2Z2GjAPOKFjShQRyU0tCuAQwpvA6EYeGtu+5YiI5A8dCSciEokCWEQkEgWwiEgkCmARkUgUwCIikSiARUQiUQCLiESiABYRiUQBLCISiQJYRCQSBbCISCQKYBGRSBTAIiKRKIBFRCJRAIuIRKIAFhGJRAEsIhKJAlhEJBIFsIhIJApgEZFIFMAiIpEogEVEIlEAi4hEogAWEYlEASwiEokCWEQkEgWwiEgkCmARkUgUwCIikSiARUQiaVEAm9lFZvaOmc0ws/vNrMTMtjGzKWY2x8weMLPiji5WRCSXNBvAZjYUOB8YHUIYARQCJwLXADeEELYHVgCndWShIiK5pqVDEEVANzMrAkqBxcAhwMTk8buBr7Z/eSIiuavZAA4hLASuAz7Gg3clMA34PIRQlWy2ABja2PPN7Awzm2pmU5cvX94+VYuI5ICWDEH0BcYB2wBDgO7AkS19gRDC7SGE0SGE0QMGDGh1oSIiuaYlQxCHAh+GEJaHECqBPwP7AX2SIQmAYcDCDqpRRCQntSSAPwb2NrNSMzNgLPAu8DxwfLLNeOCxjilRRCQ3tWQMeAr+ZdsbwNvJc24HLgW+b2ZzgC2AOzuwThGRnFPU/CYQQvgp8NMGd88FvtzuFYmI5AkdCSciEokCWEQkEgWwiEgkCmARkUgUwCIikSiARUQiUQCLiESiABYRiUQBLCISiQJYRCQSBbCISCQKYBGRSBTAIiKRKIBFRCJRAIuIRKIAFhGJRAEsIhKJAlhEJBIFsIhIJApgEZFIFMAiIpEogEVEIlEAi4hEogAWEYlEASwiEokCWEQkEgWwiEgkCmARkUgUwCIikSiARUQiUQCLiESiABYRiUQBnGf23hsuuCB2FSICYCGEznsxs+XAvE57QdkcW4cQBsQuQiSfdGoAi4hIHQ1BiIhEogAWEYlEASwiEokCWEQkEgWwiEgkCmARkUgUwCIikSiARUQiUQCLiESiABYRiUQBLCISiQJYRCQSBbCISCQKYBGRSBTAIiKRKIBFRCJRAIuIRKIAFhGJRAEsIhKJAlhEJBIFsIhIJApgEZFI/j/fTeIW/P11/gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x360 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#%% Sinkhorn\n",
    "\n",
    "lambd = 2e-3\n",
    "Gs = ot.sinkhorn(a, b, M, lambd, verbose=True)\n",
    "\n",
    "pl.figure(4, figsize=(5, 5))\n",
    "ot.plot.plot1D_mat(a, b, Gs, 'OT matrix Sinkhorn')\n",
    "\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solve Smooth OT\n",
    "--------------\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFgCAYAAACFYaNMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XmYXFW19/Hv6jnzQEISMtAgk0BQQjBhFAgCAt74igOIEq7KPIuCoFfhekFRZLgIKIoIMlwkCGhkEBAEgUQSghCmEAKZSEgCmRPS037/WLu6O51u0vOu4fd5nvNUzqlTp1ZVp1at2mefvS2EgIiIdL+i1AGIiBQqJWARkUSUgEVEElECFhFJRAlYRCQRJWARkUSUgEUiMzvAzN5IHUdnMrMTzeyfqeOQ5ikBS6eJH/aXzWy9mS0xsxvNrH+871dmtjYuVWZW3Wj9oW6ILZjZDh+1Twjh6RDCzu08/v5m9qyZrTKzD8zsGTPbu33Rto+ZVcbXWdLB43Ta39HMDjKzhY3Wy8zsT/H96Wtml5jZ7R2JN5cpAUunMLPzgSuA7wL9gPHAtsCjZlYWQjg1hNA7hNAbuBy4O7MeQvhsushdR5KWmfUFpgDXAQOB4cClwMbOia77dOXf0czKgT8B/YHDQgirOxBnh75ksoUSsHRYTECXAmeFEB4OIVSHEN4BvgxUAl9rxzEPMrOFZnaBmS01s8Vm9nkzO9LMZscq8+JG+3/KzJ4zs5Vx31+aWVm876m4279jpfaVRse/0MyWALc0rtbM7GPxOcbE9W3MbJmZHdRMuDsBhBDuCiHUhhA2hBD+FkJ4KT72xFjxXR3jm2tm+8btC+Lrm9TotfQzs9vi880zsx+YWVG8ryiuz4uPu83M+sWHZl7nyvg692l0zCvNbIWZvW1mzSbKrvg7Njp2T+AvQAlwVAhhXTuOEczsDDN7E3gzbtvFzB6Nf6s3zOzLjfbfysz+Ymarzex5M/sfy7LmGCVg6Qz7AhV4dVMvhLAWeBD4TDuPOzQedzjwQ+A3eBLYCzgA+C8z2y7uWwucBwwC9gEmAKfHOA6M+3wiVmp3Nzr+QLzCO7lJ7G8BFwK3x+RxC3BrCOHJZuKcDdSa2a1m9lkzG9DMPuOAl4CtgDuB/wP2BnaIr+mXZtY77nsdXn1uD3waOAH4z3jfiXE5ON7fG/hlvC/zOvvH1/lco+d+I743PwNuNjNrJsau+juWAw8BHwITQwgb2nkcgM/jr2dXM+sFPIq/n1sDxwI3mNmucd/rgXX433lSXLKKErB0hkHA8hBCTTP3LY73t0c1cFkIoRpPWIOAa0MIa0IIrwCvAp8ACCHMCCFMDSHUxKrt13jy+ih1wI9CCBubSwohhN8Ac4BpwDDg+80dJP6U3h8I+JfEMjP7s5kNabTb2yGEW0IItcDdwEjgv+Nz/w2oAnYws2I8kVwUX+c7wC+Ar8fjHA9cFUKYGxPjRcCxW/hJPi+E8Jv43LfG1zKkmf266u/YB/9SvDWE0NFmmZ+EED6If6+jgXfi+1oTQpgJ3At8Kb6Px+B/3/UhhFfx155VlIClMywHBrWQBIbF+9vj/Zg0ADIJ8r1G92/AK0DMbCczmxJPGq3G2ye3lDCWhRA+3MI+vwF2B677qOQRQngthHBiCGFE3H8b4JpGuzSNmxBCc69lEFAKzGt03zz8VwDxuE3vK6H5hJqxpFGc6+M/ezezX1f9HZfjXyq3mtnh7TxGxoJG/94WGBebdVaa2Ur8C2ooMBh/Xxa08NisoAQsneE5/ITTFxpvjD+pPws83g0x3Ai8DuwYQugLXAw09zO7sY8cCjDGfw1wM3CJmQ1sTSAhhNeB3+OJuK2W45X/to22jQIWxX+/28x9NXiC7+jQhl32dwwh/Ak4CZhsZgd3IMbGr3EB8I8QQv9GS+8QwmnAMvx9GdFo/5EdeN4uoQQsHRZCWIWfvLnOzI4ws1IzqwT+CCwE/tANYfQBVgNrzWwX4LQm97+Ht5m2xbXA9BDCt4C/Ar9qbqd4Iuh8MxsR10cCxwFT2/h8xIr/j8BlZtbHzLYFvg1kumrdBZxnZtvFxJjpiVCDJ5062v46M8/dpX/HEMJdwJnAA2a2X6O7isysotFS3spDTgF2MrOvx1hLzWxvM/t4fB//hH9x9oz/J07oSPxdQQlYOkUI4Wd41Xklngin4RXKhE5o92uN7wBfBdbgzQZ3N7n/Evwn8MrGZ8pbYmYTgSNoSOTfBsaY2fHN7L4GPzE0zczW4Yl3FnB+O14HwFn4yaO5wD/xk0y/i/f9Dk+ETwFv4ye2zoL65oXLgGfi6xzf1ifu6r9jCOFW/H35q5l9Km4+Dm+CySxvtfJYa4DD8OaNd/Gmlivwk37gyb5f3P4H/Mur/jWY2Sst/D27jWlAdhEpBGZ2BTA0hJA1vSFUAYtIXopNQ3uY+xTwTeC+1HE1lhdXk4iINKMP3uywDX4O4BfAA0kjakJNECIiiagJQkQkETVBCACDBg0KlZWVqcMQySkzZsxYHkIY3N7HKwELAJWVlUyfPj11GCI5xczmbXmvlqkJQkQkESVgkUKyahU8+yzMmwc6AZ+cErBIvgsB7rkHdt8d+veH/faDykrYemv43vdgzZrUERYsJWCRfLZkCUyYAF/+MhQXw2WXwQMPwA03wMEHwxVXwM47wyOPpI60IOkknEi+WrjQk+/ChXD99XDKKZ6EM047DaZNg5NOgs99Du66C445Jl28BUgVsEg+WrQIDjjAK+C//Q1OP33T5Jsxbhw8/TTsvbdXyffc0/2xFjAlYJF8U13tyXT5cnj8cW/z/Sj9+nkTxD77wKRJMGtW98QpSsAieeeCC7ynw803w9ixrXtM794webIn42OOgdXtnrBY2kAJWCSfTJkC11wDZ5/tVXBbDB0Kd98Nb70FZ5zRNfHJJpSARfLF2rXe1rv77vDzn7fvGAceCBdfDLffDo8+2rnxyWaUgEXyxSWXwIIF8OtfQ1lZ+49z8cWw446ezD/c0pyl0hFKwCL54KWXvOnhpJNg3307dqyKCu8nPGcO/OQnnROfNEsJWCQfXHAB9O0LP/1p5xzv0EPhK1/xpox33+2cY8pmlIBFct0TT3g3sosvhoEDO++4l18ONTXw3//deceUTSgBi+SyEHw8hxEj4MwzO/fY22/vV8/99rcwe3bnHlsAJWCR3HbfffCvf8Gll3rbbWf7wQ/8uP/1X51/bFECFslZIcD//I/3WDjhhK55jiFD4Nxz/RLl117rmucoYErAIrnqoYdg5ky46CIo6cJxtc49F3r06LwTfFJPCVgkF2Wq31Gj4Gtf69rnGjQITj4Z7rgD5s7t2ucqMErAIrnoH/+A557z7melpV3/fN/5jo+m9rOfdf1zFRAlYJFc9LOf+YwW3/hG9zzf8OE+Utrvfw9Ll3bPcxYAJWCRXPPqq97+e+aZ3jbbXb79bdi40a+Sk06hBCySa665xruGnXpq9z7vLrvAUUd5At6woXufO08pAYvkkqVL4bbbvNvZ4MHd//znnw/LlvkJOekwJWCRXPKrX3kzwHnnpXn+gw6CT34SrrpK09p3AiVgkVxRVQU33gif/aw3B6Rg5sn/tdd8uiPpECVgkVwxebJPsnnWWWnj+MpXvPnjuuvSxpEHlIBFcsV11/llx4cfnjaO8nK/MOMvf4G3304bS45TAhbJBdOnw9Sp3vWsKAs+tqed5nGoS1qHZMFfUkS26LrrfObiE09MHYkbPtxnT/7tb2H9+tTR5CwlYJFst3y5z1Z8wgk+60W2OPNMWLkS7rordSQ5SwlYJNvdfLN3PTv99NSRbGr//WH0aLj+enVJayclYJFsVlvrXc8OOgh22y11NJsy8y+FmTO9fVraTAlYJJs9+CDMmwdnnJE6kuZ97WveLPLLX6aOJCcpAYtks+uvh222gYkTU0fSvN69fZS0e+6B995LHU3OUQIWyVZvvumzHZ9ySveM+dtep58O1dXeI0LaRAlYJFvdeKNPNXTSSakj+Wi77AITJsCvf+3T2EurKQGLZKP16+GWW7yv7bBhqaPZsjPOgAULYMqU1JHkFCVgkWx0553exzZbT7419bnPwciR3mYtraYELJJtQvBENnq097XNBSUlPkD8Y4/B66+njiZnKAGLZJtnnoEXX/QrzcxSR9N63/oWlJWpS1obKAGLZJvrroP+/eH441NH0jZbbw3HHgu33gqrV6eOJicoAYtkk0WL4N574ZvfhF69UkfTdmedBWvX+uzJskVKwCLZ5Fe/grq67Bv3obXGjoVx47wZoq4udTRZTwlYJFt8+KH3pT36aNh++9TRtN9ZZ/lFJA8/nDqSrKcELJIt7rjDZxw+99zUkXTMl77kfZevvjp1JFlPCVgkG4TgCWuPPeDgg1NH0zFlZd6D47HH4OWXU0eT1ZSARbLBY4/BK6/4jMO51PWsJaecAj16wLXXpo4kqykBi2SDq6+GIUPguONSR9I5ttrKZ/C4/XZYujR1NFlLCVgktZdfhoce8suOy8tTR9N5zjsPqqo0ff1HUAIWSe1nP/M+v7ky7kNr7byzj2N8/fXeN1g2owQsktK8eT6p5UknwcCBqaPpfBdeCCtWwG9+kzqSrKQELJLSVVf5Sbdvfzt1JF1j/Hj49Kf9dVZVpY4m6ygBi6SydKlXhscf70M55qsLL4SFC+EPf0gdSdZRAhZJ5ec/9+nmL7oodSRd64gjYK+94LLLfOoiqacELJLC0qVwww3w1a/6yap8ZgaXXAJvv60quAklYJEUrrzSx374wQ9SR9I9jjpKVXAzlIBFutuSJd4167jj8r/6zchUwXPnaqjKRpSARbrbpZd6j4BLLkkdSfc66ijvFXHJJT7pqCgBi3Sr2bO958Mpp8AOO6SOpnuZ+UUn776rMSIiJWCR7nTxxT5IzQ9/mDqSNA44wGdQ/ulP4f33U0eTnBKwSHd56imfbug73/H50wrVT37ilyb/6EepI0lOCVikO9TU+Bi5224L3/1u6mjS2m03n3Lpxht99ucCpgQs0h2uv95HPbv6aujZM3U06f34xz5k5RlnFPTccUrAIl3t3Xe9zfeww+Dzn08dTXbo3x+uuAKefdansS9QSsAiXSkE7/FQVeUzBefDbBedZdIkPyl33nmwaFHqaJJQAhbpSrffDlOmwOWXw447po4muxQVwe9+519OJ5/sX1YFRglYpKssXAhnnw377ee3srkddvBeEQ8+6Mm4wCgBi3SF6mo49ljv/XDLLVBcnDqi7HXWWT4T9Fln+cSkBUQJWKQrfP/78MwzftWbmh4+WlER3HEH9O0LX/pSQU1fpAQs0tnuucfH+j31VK+CZcuGDYM774TXX4dvfKNguqYpAYt0pmefha9/Hfbd1/v8SusdcoiPFXHPPf4LogCUpA5AJG+88YbPAjxyJDzwAFRUpI4o95x/Prz1lo8VMWoUnHZa6oi6lBKwSGd4/XWv4Mz8jP6gQakjyk1mcN113oPk9NP95OXJJ6eOqsuoCUKko2bN8rP4dXXw5JM66dZRJSUwebKPH3zKKX4BS55SAhbpiIce8vZeM3jiCdh119QR5Yfych85buJE75523nlQW5s6qk6nBCzSHrW1fnXb0Uf7xQT/+hd8/OOpo8ovmSR87rlwzTVw5JE+nVMeUQIWaau5c73J4fvf936rTz0FI0akjio/FRd7b5KbbvL3efRouP/+1FF1GiVgkdZat85HNdt1Vx/H9rbb4K67oHfv1JHlv5NOghkz/Ivu//0/bx+ePTt1VB2mBCyyJatWebeo7bbzcWy/8AV49VXv76vRzbrPrrvCtGlw5ZXw9NO+fsIJ/rfIUUrAIs2pqYG//x1OPNGv0rroIthrL7+8+M471eSQSlmZ9xWePRvOOcfbiHfbDQ480Ke7X7EidYRtYqEAh4CTzY0dOzZMnz49dRjp1Nb6hRTPPOO9GR55BD74wJsXvvpV7w41ZkzqKKWpZct8FLWbb4Y33/Q24wMPhAkT/HbMGOjVq8ue3sxmhBDGtvvxSsACBZCAq6q8Olq2DN57zzv6z5sHc+Z44p01C9av932HDoVDD/W2xsMP79IPsHSSEOD55+G+++Cvf/Xpn8CbiHbayXuo7LgjVFb6lYrDhvnEqAMH+t+3nU1JSsDSKcYOHhymT5zY9U/U0v+3zPbG94ew+VJX57e1tQ1LTY0v1dWwcaMvGzb4snYtrFnj/27OiBGw886w++5eLY0b5x9Yte3mtvff93E5Zs70E6ZvvOFftlVVm+9bUgJ9+vjSsyf06OFd4MrLvcmjtNT3KS5uWLbZBq66SglYOsfYsrIwvbumSm8puWW2N77fbNOlqMhvG38Yiov9Q1Ja2vDB6dHDlz59vBmhf39fBg+GIUM88Y4Y4ftKYair818/8+f77dKl/qtoxQr/kl6zxn8FbdjQ8EVeVeVf7DU1DV/4dXU+u/Ujj3Q4AWssCHF77AH53AQhUlTkTQ/DhqWOpJ56QYiIJKIELCKSiNqABQAzWwO8kTqOLRgELE8dxBbkQoyQG3HmQow7hxD6tPfBagOWjDc6cjKhO5jZdMXYOXIhzlyJsSOPVxOEiEgiSsAiIokoAUvGTakDaAXF2HlyIc68j1En4UREElEFLCKSiBKwiEgiSsAFzsyOMLM3zGyOmX0vdTwZZjbSzJ4ws1fN7BUzOyduH2hmj5rZm/F2QBbEWmxmM81sSlzfzsymxff0bjMrSxxffzObbGavm9lrZrZPtr2PZnZe/DvPMrO7zKwiG95HM/udmS01s1mNtjX73pn73xjvS2a2xfFLlYALmJkVA9cDnwV2BY4zs2yZ1rcGOD+EsCswHjgjxvY94PEQwo7A43E9tXOA1xqtXwFcHULYAVgBfDNJVA2uBR4OIewCfAKPNWveRzMbDpwNjA0h7A4UA8eSHe/j74Ejmmxr6b37LLBjXE4Gbtzi0UMIWgp0AfYBHmm0fhFwUeq4Woj1AeAz+NV6w+K2YfgFJCnjGhE/hIcAUwDDr94qae49ThBfP+Bt4gn3Rtuz5n0EhgMLgIH4xWFTgMOz5X0EKoFZW3rvgF8DxzW3X0uLKuDClvmPn7EwbssqZlYJ7AlMA4aEEBbHu5YAQxKFlXENcAFQF9e3AlaGEGrieur3dDtgGXBLbCb5rZn1IovexxDCIuBKYD6wGFgFzCC73sfGWnrv2vx5UgKWrGZmvYF7gXNDCKsb3xe8zEjWj9LMjgaWhhBmpIqhFUqAMcCNIYQ9gXU0aW7IgvdxADAR/7LYBujF5j/7s1JH3zsl4MK2CBjZaH1E3JYVzKwUT753hBD+FDe/Z2bD4v3DgKWp4gP2A/7DzN4B/g9vhrgW6G9mmXFWUr+nC4GFIYRpcX0ynpCz6X08FHg7hLAshFAN/Al/b7PpfWyspfeuzZ8nJeDC9jywYzzbXIaf+Phz4pgAP6MM3Ay8FkK4qtFdfwYmxX9PwtuGkwghXBRCGBFCqMTfu7+HEI4HngC+GHdLHeMSYIGZ7Rw3TQBeJYveR7zpYbyZ9Yx/90yMWfM+NtHSe/dn4ITYG2I8sKpRU0XzUjW8a8mOBTgSmA28BXw/dTyN4tof/2n3EvBiXI7E21gfB94EHgMGpo41xnsQMCX+e3vgX8Ac4B6gPHFsnwSmx/fyfmBAtr2PwKXA68As4A9AeTa8j8BdeLt0Nf5r4pstvXf4Cdjr42fpZbxXx0cev0OXIpvZEfhPrmLgtyGEn7b7YCIiBabdCTj2IZ2Ndw1aiP+cPS6E8GrnhScikr86MiD7p4A5IYS5AGb2f/iZzBYT8KBBg0JlZWUHnlI6auVKn/x15MhNt8+YMWN5CGFwZv0zRV/SKE0ijTxad08L03m3X0cScHN93sY13cnMTsavCmHUqFFM18y7SV1wAVx33eYTIJvZvDQRiRSuLu8FEUK4KYQwNoQwdvDgwVt+gHSp6mooLU0dhYhAxxJwVvchleZVVUFZ0qFhRCSjIwk4a/uQSsvWr4devVJHISLQgTbgEEKNmZ0JPIJ3Q/tdCOGVTotMusTq1dC7d+ooRAQ6OC19COFB4MFOikW6wfvvw1ZbpY4ii1g8sd2B/vAi7aVLkQvMokUwdGjqKEQElIALSm0tvPMOfOxjqSNJyKzJUhSXptubLCJdQAm4gMyaBTU1sPvuqSMREehgG7Dkln/8w28PPDBtHNnEiry6DbWb3eG3IY6z3lIVrLZj6QBVwAXk9tth1103vwxZRNJQBVwgpk6F55/3y5ALmm1ac1hJCx+B4mK/rYsVbqyU69czlTGZ1ea3q0KWj6IKuABUVcFpp8GQIXDCCamjEZEMVcB5LgT4/vfhxRfh/vuhb9/UEWWXzYZjjRWyZSrg4sz2Jm3ATR5ndXWbHq9JRbxZhazKWFAFnNdC8NHPrrzSK+CJE1NHJCKNqQLOU6tWwdlnw223wRlnwP/+b+qIskTTNtq6FirRJhWvlcaPSpM25Pq24cx+tbWbHDfE9cz2UJupgFUZiyrgvHT//d7b4fbb4Yc/9BNvRfpLi2QdVcB55Lnn4Cc/gb/8BfbYwxPx3nunjirLZCpMy/T/3bQDsBXVNbtffeVbnLn1xmHLDK7cdEyJuG41Nb45thHXr1dV+3p9ZRzjyKyrMi4IqotyXG0t3Hcf7Lcf7Lsv/POfcPnlPuOFkq9IdlMFnKNeew3uvhvuuAPmzIHKSm/n/c//1HCTrVJfUTbtzxtrkkzbbX1lm9mvaNPtJV4Jh7JNK+FQXrLJ4TO9JKjx4xbFCjizHjZW+frGjb6eqZSr422mQlZFnFeUgHPI7NmedP/4Rx/XwcwvK77sMvjCF6ClawpEJDvpI5vF1q6Fp5+Gxx6DRx+Fl1/2pLv//n5i7ZhjYNiw1FHmuM0qydg7IVae9VfEZdpsM70b4t4W4vxOTSrf2l7lvl7iFXMojm3C8fmsum6T4xVtiMffUBXXYyX8od/WV8ZVfn+mN4Uq49ymBJxFqqth2jRPuI8/7pcP19T4HG777QfXXANf/CIMH546UhHpDErACS1d6gk3szz3HKxb58XUXnvB+efDoYd68u3RI3W0BaJJ23CoadpvuEmlmekNkWkTbtIvuK7MK+CaHr5fbXmT+0t9vXijH6/4Qz9+6XqvwIvXxop4zYd+u95vw4YNfvth823GqohzgxJwN9m40S8HnjrVk+3UqfD2235fcTGMHg2TJsGECXDQQTBwYNJwRaQbKAF3gbVr4aWXPOHOnOm3L73kg+KANyGMH++XB48b59WuZirOMk0qx1ATey2E4mbvz1z5ZrFttrg07hcr4tpyr4TrSmNFXOHbq/rGtuNMd+M6f1xxlX80S9ZXAFC+xv+DlK30CrdklVe+RWvW++1avw3rY2WcaSuuadqLQhVxNlEC7qD33ts00c6cCW++2fD/fOBA2HNPOOccT7bjxsGIEWljFpHsoATcShs2wKuvek+ExsuSJQ37VFZ6sj3+ePjkJ/3fI0ZoSrG8EDYd26F+DInMbWZ7rDgzVzhZfRtyT398rIjrSjIVsa9Xx19AtT3i88QC2uLDizf4R7Vstd+Wr/BeFhUr/IEVy73iLVkRK+LV6/w4mbbiDbHtWP2Ks4oScBO1tTB37uaJds6chvMvFRU+1sLhhzck2k98Avr3Txu7iOSWgk7AS5dunmhfeQXWexGBmc8gPHo0HHus344eDTvs0NA9VApMplKMk8iF0GR0s1hZFmWupKv2tuOSWHlaLHWLarxttyb2isj0E66JlXBN/9gLolfs1RDHqFhb45Xz6rV+5V3ZB/4fsWKZV9g9l/pxeyzt4/e/75Vv8cq1/jxr4m2mIo79m1URp1EQCXj9ek+sTZPt0qUN+wwe7Mn1pJMaEu1uu+nkmIh0nbxLwB984CfCXnih4Xb27IYv9IoKn5b9qKMaEu3o0T5dj0ibtdA2XJfpFVGTqXxjJRwrzqIP/Zvd6nzgjvUb40cxeIW7IY4xQayAt+rvlevgnt62W2L+fKuqvIP4klVe8S5ZGtuEF3sbca/FZfE2bl/qU6IUf+DHszWxrbhp7wn1J+4WOZuAQ4B339082c6f37DPyJEwZsymzQcf+5iaD0QkO+RMAq6r82aEp57y8RGeegoWL264f6edYJ99fPaHPff0ZdCgdPFKgWmhbbjpjBiWGVMiXsFWsTH2613rFWpp7PdbssGrhLXVvr4sPk3fMn/czgPeA2BU+Qd+RxwTZP52fgXPyyu3AWDOksEArJrvlXLvhV4p937X24x7LPG24JLlXhEXrfbbunWxX7Eq4i6VtQm4uhpmzGhIts88AytW+H3Dh/vVYuPHe4X7iU9Anz5JwxURabOsSsAh+IDit93mQy6uXu3bd9rJh1s84AAffrGyUn1rJcs1rYirYuWYaSuObcKZirh0nbfBFq+OvRdWxSvfVnsb7trVXsHOXuel7rpRvr10iB9vfO85ABzQczYAH/b3j/br2/j+T2+/EwDT3t0WgHffiZXwfH+ePgu90u652LeXLFsDqCLualtMwGY2ErgNGIKPwndTCOFaMxsI3A1UAu8AXw4hrGhPEHPnwh/+4Il37lzveXDMMfC5z/nQi0OHtueoIiLZrTUVcA1wfgjhBTPrA8wws0eBE4HHQwg/NbPvAd8DLmxrAL/4BXznO17RTpgAl1zi1a66f0leyfSWaDI2Q2bsCIvj/Vq8cq1inVeipbF3Q8UH3ma75gP/yC5937vtTN7Oe1HMH+Vtv0dv9W8AxlUsAGCP3gsBOLznXAD+PWgrAB7edjQA/1i4AwCL3/ariHq/48/TZ4H3oui52I+fqYgtVsT1vSYyo7GpH3G7bDEBhxAWA4vjv9eY2WvAcGAicFDc7VbgSdqYgG+6yZPvMcfA1Vd7rwURkULRpjZgM6sE9gSmAUNicgZYgjdRtNrbb8Opp8LYsXDnnT7ouEjBqK+IN70SzTJXzMXKMjPKWa+VXomWv+/9eHvFK95WL/Htzy3eBYBXtvf2ukNGeFvw5/q/CMCe5X78I3r6cceVPwXAtH4vAzBlm08C8OQor4iXzPXKu887mTZir4h7vOvrJcu9Ig71/YgzbcTNXFmnarhFrZ4V2cx6A/dlzgyZAAAPx0lEQVQC54YQVje+L4QQaJilpenjTjaz6WY2fdmyZfXbR4yAgw/2/rsPPti+4EVEcpmFVnw7mVkpMAV4JIRwVdz2BnBQCGGxmQ0Dngwh7PxRxxk7dmyYPn16/fqaNfCZz8Dzz/vANpMmwcSJfrWadC8zmxFCGJtZ/0zRl1S2pJDp3mNx9uVS/5Fq8Sei9YxTo/TzCrVmsN+uH+YV6pqRsf/wtt4W27tyFQCfHuG9JA7t9woAu5X5dfg949O9W+vHf3b9jgA8uvzjAMya5/2Jy97x4/ee5/v3WeSVesUSr4CLVsSKeG1mFLbMWBNVeTMW8aN193R636stVsDmc63cDLyWSb7Rn4FJ8d+TgAfa+uR9+sDDD8OFF/rYDMce6z0eTj7Z+/3W1W35GCIiuWqLFbCZ7Q88DbxMZqIsuBhvB/4jMAqYh3dD++CjjtW0Am6sthaefBJuvRXuvdcH0BkwwLuhHXig9wEeMwZKS9vy8qS1VAFnqRYq4qJyr0gzFXGIFXH14Di2xDb+M3LNCH/culH+0e0VK+K9h3oviX37eWW8Y7kPbF0cP+ILqr23xLNrvE34uSXbAbB8vveW6DXf4+gzPx53sfcPLlvqFbCtjBXxunWbz1uXoxVxV1TArekF8U+gpSee0FmBFBd7N7QJE+CGG+CBB+CJJ/xKuL/8xffp2dMvN85ckPGpT6m7mojkrqy6Ei6jd2+fVeL44319yRK/Qi4zDsSll/qXpxnsvLOP+zBmTMMYEJrQUvJGC1fU1WZ6S8T+w0WxX27ZKq88S5d5RdxrUayI53vFvHbuAACeGtE/3nqF+7GhfoJ8t37esWl4+UoARlb4j9q6IV5JzyrzSnd+L/+QVffx427sF9uI+3pK6RFHYSv+oLx+xDVijKjvcL2sTMBNDR0KX/yiLwArV8Kzz/rJu5kzPTnfdVfD/ttu25CQx4zxWSu22UaXL4tIdsmJBNxU//5w5JG+ZCxf7sk4s7zwAtx/f8OX6oABDUNS7rGH3+6+uwbxkRzTyorY4tgNJSu9x2jfpV4J91rotxvmeRvxuiHehjxvmI8R8caw4QD0HOxV69Z9/cq3fmXeq6FnqVfAW23l29+PRc2aUq+Aa3qUxFtvG+zZq5TS9/2+ohXxKrq1Teary/G24Y7IyQTcnEGDvEvbZz7TsG3NGvj3v315+WWfGv6223x7RmXlpgOzjx7tg//oZJ/khBYScmawHIuD51j8+V8cT471iQO391zotxsHeSLeMNhTwvqt+wGwaJBf+PHOwDjVUi+/0KK4JF5CXeTPX9PX7/+wLjObqDdZ1JaV06NHnDapLHMC0T9ctiqmnw89uVOAA/3kTQJuTp8+3oti//0btoUA8+ZtPj3RQw/VT2hLWRnsssvmiVkzHItIZ8rrBNwcM696Kyt9tLWMjRvhjTe8Ss4k5X/8A+64o2Gf/v292aJpU0bfvt39KkRasKWKOJ4As3jpcNEqb6LoudQH4ano500UfQZ4RbxxK69WNwz0KnbjAE8Z1bHpLlTEqZcys8zEp49zjrKxvxEsk2b8mGWlvnNJnJqmKE4wGtbFpglr0m0tTu9EXW1r34WcUXAJuCXl5Z5Q99hj0+0rVsCsWZtWy3fc0TBWMfgsyZnp6TO3Q4eqWhaRj6YEvAUDBni/4wMOaNgWAixY4NXyiy/68sILMHlywz5bb92QkD/5SdhrL0/USsrSrbY0VVKTk3YWT5CUve/VaumSON19X1+v7uelbVUc8H1jH2/rrfHdqC3b9D+41UJdHGiruncsk0O8rDrGVlzkj7GieGFubD8mxkas3utbgvOobVgJuB3MYNQoX44+umH7qlWelGfO9KQ8cyZcdVX95AcMHAjjxvkyfrxfSDJgQJrXICLpKQF3on79Nq+Wq6rg1Vdh+nSYNg2mTvXxLzJf3jvt5Ml43DjYd19vAilq9Rh1Im3UtCLOFJOZijhzJjoOpmNrM4PteNVa0cMr4fLefturt1fENb29Hbemp1e5NRX+n7iuxAiZwrY2DsFZ7BVvbUUcaKgmdlOLIVpmEJjiTCXsvSMyvSQy0zrlQ7c1JeAuVlbW0AzxrW/5ttWrPSFPnepJ+eGHvXsceHe6gw+GQw/1y7K3317NFiL5Sgk4gb594ZBDfIGGrnFPPQWPP+7LPff4fZWVDWNkHHGEmiykizSZMqm+ytys94T3J7Y1XvFmBgUqr/DbsngbenjFXNujlLryTBeJTRVVx0o3Vhgh9o6wsthPuOkD6jbdvz70HK6ElYCzQOOucSec4P+P3nijIRnfey/cfLNfHHLYYfDlL/u4yf36pY5cRDpCCTgLmfmFILvsAmec4cXI8897Iv7jH+Gvf/WmjSOO8GR8zDEaxF46WUu9J+IUSpkr7OrbikuaDBwfq9jSsjJCRewGUeIVboj9f+sbfTOFbdMKN+5vdfGy1DiBaX11HmPM1MP1lXDT15DFdLonBxQX+4m6n/8c3nkHnnvOE/OMGfC1r/ngQz/+Mbz/fupIRaQtWjUlUWf5qAHZpe3q6nzM5F/8wi+l7tEDvvEN+O53PSm3hQZkl3ZpOmB87NNLrHKtpKR+EHkyVXLJpuv1x8g8trZJP99YGYdM1R2vjAsbNx3Wsn4MiYxO7i+cZEoiyV5FRX5y7sEHG6Z0uukm2HVXuOaa+l9qIpKllIDzxO67w+9+B2++CQcdBOed5/2KX3kldWSS10Lwpa4W6moJtXGpqiJUVVG34UPq1m3wZc1aX1av8WXVaupWrSasWePL2nW+bNjgy8aNvlRXE6qrvaKorY1T3TdMGGlmvhRtujTaIWv7cioB55ltt4UpU+DOO2HuXE/CzzyTOioRaY56QeQhMzjuOB+G89BDvevan//szRUiXappe2uore9BUb+pSXtx/Xpxk3qwfmwI2/TYmfVM23BLbbzW5Hgh+9rkVAHnsZEj/eKO7bf36ZwWLUodkYg0pgSc54YMgfvu84GlTjopJ7pGSr7JtBO31F5cU02oqaauatMl044cqqp9qa7xng6xLTjU1hFq67xXRF1u/sdWAi4AO+wAl1/uXdWeeCJ1NCKSoQRcIE491YfDvOGG1JGIRC1Uxi1VyE0r5fr7M0td2GRp2lsiGykBF4iKCpg0CR54AOLMLyKSmBJwATnsMJ94dOrU1JGItEKmMm6j+n7AVgRWtHm/4CyiBFxA9tnHb//1r7RxiIhTP+AC0q+f94p4663UkYi0QTN9iz9y97qPXs8mqoALTGUlzJ+fOgoRASXggjNokIatFMkWSsAFZsAAWLEidRQiAm1IwGZWbGYzzWxKXN/OzKaZ2Rwzu9vMyrouTOksPXvChg2poxARaFsFfA7wWqP1K4CrQwg7ACuAb3ZmYNI1KiqUgEWyRasSsJmNAI4CfhvXDTgEmBx3uRX4fFcEKJ2rtBSqq1NHISLQ+gr4GuAC6qfPYytgZQghMwfIQmB4cw80s5PNbLqZTV+2bFmHgpWOKynxizFEJL0tJmAzOxpYGkKY0Z4nCCHcFEIYG0IYO3jw4PYcQjpRUVH9MKoiklhrLsTYD/gPMzsSqAD6AtcC/c2sJFbBIwCNNpsDlIBFsscWK+AQwkUhhBEhhErgWODvIYTjgSeAL8bdJgEPdFmU0mnMNCawSLboSD/gC4Fvm9kcvE345s4JSbqSErBI9mjTWBAhhCeBJ+O/5wKf6vyQpCtl6eSwIgVJV8KJiCSiBFxgVAGLZA8lYBGRRJSARUQSUQIWEUlECVhEJBElYBGRRJSARUQSUQIWEUlECVhEJBElYBGRRJSARUQSUQIWEUlECVhEJBElYBGRRJSARUQSUQIWEUlECVhEJBElYBGRRJSARUQSUQIWEUlECVhEJBElYBGRRJSARUQSUQIWEUlECVhEJBElYBGRRJSARUQSUQIWEUlECVhEJBElYBGRRFqVgM2sv5lNNrPXzew1M9vHzAaa2aNm9ma8HdDVwYqI5JPWVsDXAg+HEHYBPgG8BnwPeDyEsCPweFwXEZFW2mICNrN+wIHAzQAhhKoQwkpgInBr3O1W4PNdFaSISD5qTQW8HbAMuMXMZprZb82sFzAkhLA47rMEGNLcg83sZDObbmbTly1b1jlRi4jkgdYk4BJgDHBjCGFPYB1NmhtCCAEIzT04hHBTCGFsCGHs4MGDOxqviEjeaE0CXggsDCFMi+uT8YT8npkNA4i3S7smRBGR/LTFBBxCWAIsMLOd46YJwKvAn4FJcdsk4IEuiVBEJE+VtHK/s4A7zKwMmAv8J568/2hm3wTmAV/umhBFRPJTqxJwCOFFYGwzd03o3HBERAqHroQTEUlECVhEJBElYBGRRJSARUQSUQIWEUlECVhEJBElYBGRRJSARUQSUQIWEUlECVhEJBElYBGRRJSARUQSUQIWEUlECVhEJBElYBGRRJSARUQSUQIWEUlECVhEJBElYBGRRJSARUQSUQIWEUlECVhEJBElYBGRRJSARUQSUQIWEUlECVhEJBElYBGRRJSARUQSUQIWEUlECVhEJJFWJWAzO8/MXjGzWWZ2l5lVmNl2ZjbNzOaY2d1mVtbVwYqI5JMtJmAzGw6cDYwNIewOFAPHAlcAV4cQdgBWAN/sykBFRPJNa5sgSoAeZlYC9AQWA4cAk+P9twKf7/zwRETy1xYTcAhhEXAlMB9PvKuAGcDKEEJN3G0hMLy5x5vZyWY23cymL1u2rHOiFhHJA61pghgATAS2A7YBegFHtPYJQgg3hRDGhhDGDh48uN2Biojkm9Y0QRwKvB1CWBZCqAb+BOwH9I9NEgAjgEVdFKOISF5qTQKeD4w3s55mZsAE4FXgCeCLcZ9JwANdE6KISH5qTRvwNPxk2wvAy/ExNwEXAt82sznAVsDNXRiniEjeKdnyLhBC+BHwoyab5wKf6vSIREQKhK6EExFJRAlYRCQRJWARkUSUgEVEElECFhFJRAlYRCQRJWARkUSUgEVEElECFhFJRAlYRCQRJWARkUSUgEVEElECFhFJRAlYRCQRJWARkUSUgEWkMJj5kkWUgEVEEmnVjBgiIjkvhNQRbEYVsIhIIkrAIiKJKAGLiCSiBCwikogSsIhIIkrAIiKJKAGLiCSiBCwikogSsIhIIkrAIlIYiop9ySJKwCIiiWgsCBEpDHW1qSPYjCpgEZFElIBFRBJRAi4w48fDOeekjkJEACx04xiZZrYMmNdtTyhtsW0IYXDqIEQKSbcmYBERaaAmCBGRRJSARUQSUQIWEUlECVhEJBElYBGRRJSARUQSUQIWEUlECVhEJBElYBGRRJSARUQSUQIWEUlECVhEJBElYBGRRJSARUQSUQIWEUlECVhEJBElYBGRRJSARUQSUQIWEUlECVhEJBElYBGRRJSARUQS+f/Js7GEh8ojHQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x360 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFgCAYAAACFYaNMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8lNXZ//HPRcIqyCIIyBYpgkVwQSrg9mixSl0e/NVqtVqxteK+t1q1T6tPH61a61LqUitarUpVrEupS9XiUhUUxIVFEEEEZAn7qiHJ+f1x3TEBAmSZyZnl+3695jWZmTsz19xJvjlz7nOfYyEERESk4TWKXYCISL5SAIuIRKIAFhGJRAEsIhKJAlhEJBIFsIhIJApgkSrM7BAzmxm7jlQyszPM7D8RX/8wM1sQ6/UzmQJYUir5Y//IzDaY2WIzu9vM2iSP3WNm65JLiZltqnL7+QaoLZhZr+1tE0J4I4TQp47Pf7CZvWVmq81shZm9aWbfqlu1dWNmRcn7LKzn86Tl52hmTc1stJnNM7O1Zva+mX23PrVmMwWwpIyZXQ7cBPwcaA0MBnoAL5lZkxDCOSGEliGElsANwGMVt0MI0f8I6xNaZrYzMA4YBbQDugDXAV+lprqGk+afYyEwH/iv5Ll/CTxuZkU1rK1e/1gyjQJYUiIJoOuAC0MIL4QQNoUQPgNOAoqA0+rwnIeZ2QIzu8LMlprZIjM73syONrNZSSvz6irbH2Bmb5vZqmTbP5pZk+Sx15PNPkhaaj+o8vxXmtli4IGqH5fN7BvJawxIbu9mZsVmdlg15fYGCCGMCSGUhRA2hhD+FUL4MPneM5IW8W1JfXPM7MDk/vnJ+xtR5b20NrOHktebZ2a/NLNGyWONktvzku97yMxaJ99a8T5XJe9zSJXnvMXMVprZ3G21OtPxc6wqhLA+hHBtCOGzEEJ5CGEcMBfYfxv1VN1vy4Frk/t/YmYzkvfzopn1qPI9R5rZzOSTyF1m9pqZ/bQ+daeLAlhS5UCgGfD3qneGENYBzwHfqePzdkqetwvwK+DPeAjsDxwC/I+Z7Z5sWwZcCrQHhgBDgfOSOg5Nttknaak9VuX52+EtvJFb1P4pcCXwsJm1AB4AHgwhvFpNnbOAMjN70My+a2Ztq9lmEPAhsAvwKPA34FtAr+Q9/dHMWibbjsJbiD3x1uLpwI+Tx85ILocnj7cE/pg8VvE+2yTv8+0qrz0z2Tc3A6PNzKqpMV0/x2qZWUf8n9e07Ww2CJgDdASuN7PhwNXA94AOwBvAmOT52gNjgavw/TwTf08ZSQEsqdIeWBZCKK3msUXJ43WxCbg+hLAJD6z2wB0hhLUhhGnAdGAfgBDC5BDChBBCadJq+xMeXttTDvw6hPBVCGHjlg+GEP4MzAYmAp2Ba6p7khDCGuBgIOD/JIrN7NkkYCrMDSE8EEIoAx4DugH/m7z2v4ASoJeZFQAnA1cl7/Mz4PfAj5LnORW4NYQwJwnGq4CTd/DxfF4I4c/Jaz+YvJeO1WyXrp/jVsysMfAI/k/t4+1s+kUIYVTyc90InAP8NoQwI6nzBmDfpBV8NDAthPD35LE/AItTVXOqKYAlVZYB7bcRAp2Tx+tieRIaABUBuaTK4xvxFiBm1tvMxiUHjdbgf5g7CoziEMKXO9jmz0A/YFQIYZt9ukkgnBFC6Jpsvxtwe5VNtqybEEJ176U90BiYV+WxefinAJLn3fKxQqoP1Apfh1AIYUPyZctqtkvXz3EzSXfKX/F/OhfsYPP5W9zuAdyRdOWsAlYAhu+f3apuH3y2sYwdgaEAllR5Gz/g9L2qdyYfqb8LvNIANdwNfAzsEULYGf+YWt3H7Kq2Ox1gUv/twGjgWjNrV5NCkhbdX/Agrq1leMu/R5X7ugMLk6+/qOaxUjzg6zu9Ydp/jknXx2j8H8YJyaeb7dnyPc0Hzg4htKlyaR5CeAtvpXfd4rW6kqEUwJISIYTV+MGbUWY2zMwaJ0e2H8dbIH9tgDJaAWuAdWa2J3DuFo8vwftMa+MOYFII4afAP4F7qtvIzPY0s8vNrGtyuxtwCjChlq9H0uJ/HO/vbJV8tL4MeDjZZAxwqZntngRjxUiEUqAY71ap7fuseO2G+DneDXwTOK66bp8auAe4ysz2gq8PWJ6YPPZPoH9ysLYQOB/v589ICmBJmRDCzXir8xY8CCfirZWh2/vonkI/A34IrMW7DR7b4vFrgQeTj64n7ejJkoM9w6gM8suAAWZ2ajWbr8UPFk00s/V48E4FLq/D+wC4EFiPH3z6D37Q7v7ksfvxIHwdH0HwZbJ9RffC9cCbyfscXNsXTufPMflncjawL7DYKscPV7dPt1XfU/gwub8lXU1T8dY5IYRlwIn4gcblQF9gEslwQPMTbdbV5z2kkmlCdhHJVUlf8wLg1BDC+Nj1bEktYBHJKWZ2lJm1MbOmVB4HqHVXUENQAItIrhkCfIofzDwOOL6Ofc1ppy4IEZFI1AIWEYkkpya2kLpr3759KCoqil2GSFaZPHnyshBCh7p+vwJYACgqKmLSpEmxyxDJKmY2b8dbbZu6IEREIlEAi+ST1avhrbdg3jzQAfjoFMAiuS4EeOIJ6NcP2rSBgw6CoiLYdVf4xS9g7drYFeYtBbBILlu8GIYOhZNOgoICuP56eOYZuOsuOPxwuOkm6NMHXnwxdqV5SQfhRHLVggUevgsWwJ13wtlnewhXOPdcmDgRzjoLjjsOxoyBE06IV28eUgtYJBctXAiHHOIt4H/9C847b/PwrTBoELzxBnzrW95KfuKJhq81jymARXLNpk0epsuWwSuveJ/v9rRu7V0QQ4bAiBEwdWrD1CkKYJGcc8UVPtJh9GgYOLBm39OyJYwd62F8wgmwZk16axRAASySW8aNg9tvh4su8lZwbXTqBI89Bp9+Cuefn576ZDMKYJFcsW6d9/X26we/+13dnuPQQ+Hqq+Hhh+Gll1Jbn2xFASySK669FubPhz/9CZo0qfvzXH017LGHh/mXO1qvVOpDASySCz780LsezjoLDjywfs/VrJmPE549G37729TUJ9VSAIvkgiuugJ13hhtvTM3zHXEE/OAH3pXxxRepeU7ZigJYJNuNH+/DyK6+Gtq1S93z3nADlJbC//5v6p5TNqMAFslmIfh8Dl27wgUXpPa5e/b0s+fuuw9mzUrtcwugABbJbk89Be+8A9dd5323qfbLX/rz/s//pP65RQEskrVCgP/7Px+xcPrp6XmNjh3hkkv8FOUZM9LzGnlMASySrZ5/HqZMgauugsI0zqt1ySXQvHnqDvDJ1xTAItmoovXbvTucdlp6X6t9exg5Eh55BObMSe9r5RkFsEg2eu01ePttH37WuHH6X+9nP/PZ1G6+Of2vlUcUwCLZ6OabfUWLn/ykYV6vSxefKe0vf4GlSxvmNfOAAlgk20yf7v2/F1zgfbMN5bLL4Kuv/Cw5SQkFsEi2uf12Hxp2zjkN+7p77gnHHOMBvHFjw752jlIAi2STpUvhoYd82FmHDg3/+pdfDsXFfkBO6k0BLJJN7rnHuwEuvTTO6x92GOy7L9x6q5a1TwEFsEi2KCmBu++G737XuwNiMPPwnzHDlzuSelEAi2SLsWN9kc0LL4xbxw9+4N0fo0bFrSMHKIBFssWoUX7a8VFHxa2jaVM/MeMf/4C5c+PWkuUUwCLZYNIkmDDBh541yoA/23PP9To0JK1eMuAnKSI7NGqUr1x8xhmxK3FduvjqyffdBxs2xK4maymARTLdsmW+WvHpp/uqF5niggtg1SoYMyZ2JVlLASyS6UaP9qFn550Xu5LNHXww9O8Pd96pIWl1pAAWyWRlZT707LDDYK+9YlezOTP/pzBlivdPS60pgEUy2XPPwbx5cP75sSup3mmnebfIH/8Yu5KspAAWyWR33gm77QbDh8eupHotW/osaU88AUuWxK4m6yiARTLVJ5/4asdnn90wc/7W1XnnwaZNPiJCakUBLJKp7r7blxo666zYlWzfnnvC0KHwpz/5MvZSYwpgkUy0YQM88ICPte3cOXY1O3b++TB/PowbF7uSrKIAFslEjz7qY2wz9eDblo47Drp18z5rqTEFsEimCcGDrH9/H2ubDQoLfYL4l1+Gjz+OXU3WUACLZJo334T33/czzcxiV1NzP/0pNGmiIWm1oAAWyTSjRkGbNnDqqbErqZ1dd4WTT4YHH4Q1a2JXkxUUwCKZZOFCePJJOPNM2Gmn2NXU3oUXwrp1vnqy7JACWCST3HMPlJdn3rwPNTVwIAwa5N0Q5eWxq8l4CmCRTPHllz6W9thjoWfP2NXU3YUX+kkkL7wQu5KMpwAWyRSPPOIrDl9ySexK6ufEE33s8m23xa4k4ymARTJBCB5Ye+8Nhx8eu5r6adLER3C8/DJ89FHsajKaAlgkE7z8Mkyb5isOZ9PQs205+2xo3hzuuCN2JRlNASySCW67DTp2hFNOiV1Jauyyi6/g8fDDsHRp7GoylgJYJLaPPoLnn/fTjps2jV1N6lx6KZSUaPn67VAAi8R2880+5jdb5n2oqT59fB7jO+/0scGyFQWwSEzz5vmilmedBe3axa4m9a68ElauhD//OXYlGUkBLBLTrbf6QbfLLotdSXoMHgz/9V/+PktKYleTcRTAIrEsXeotw1NP9akcc9WVV8KCBfDXv8auJOMogEVi+d3vfLn5q66KXUl6DRsG++8P11/vSxfJ1xTAIjEsXQp33QU//KEfrMplZnDttTB3rlrBW1AAi8Rwyy0+98Mvfxm7koZxzDFqBVdDASzS0BYv9qFZp5yS+63fChWt4DlzNFVlFQpgkYZ23XU+IuDaa2NX0rCOOcZHRVx7rS86KgpgkQY1a5aPfDj7bOjVK3Y1DcvMTzr54gvNEZFQAIs0pKuv9klqfvWr2JXEccghvoLyjTfC8uWxq4lOASzSUF5/3Zcb+tnPfP20fPXb3/qpyb/+dexKolMAizSE0lKfI7dHD/j5z2NXE9dee/mSS3ff7as/5zEFsEhDuPNOn/XsttugRYvY1cT3m9/4lJXnn5/Xa8cpgEXS7YsvvM/3yCPh+ONjV5MZ2rSBm26Ct97yZezzlAJYJJ1C8BEPJSW+UnAurHaRKiNG+EG5Sy+FhQtjVxOFAlgknR5+GMaNgxtugD32iF1NZmnUCO6/3/85jRzp/6zyjAJYJF0WLICLLoKDDvJr2VqvXj4q4rnnPIzzjAJYJB02bYKTT/bRDw88AAUFsSvKXBde6CtBX3ihL0yaRxTAIulwzTXw5pt+1pu6HravUSN45BHYeWc48cS8Wr5IASySak884XP9nnOOt4Jlxzp3hkcfhY8/hp/8JG+GpimARVLprbfgRz+CAw/0Mb9Sc9/+ts8V8cQT/gkiDxTGLkAkZ8yc6asAd+sGzzwDzZrFrij7XH45fPqpzxXRvTuce27sitJKASySCh9/7C04Mz+i37597IqykxmMGuUjSM47zw9ejhwZu6q0UReESH1NnepH8cvL4dVXddCtvgoLYexYnz/47LP9BJYcpQAWqY/nn/f+XjMYPx769o1dUW5o2tRnjhs+3IenXXoplJXFrirlFMAidVFW5me3HXusn0zwzjvwzW/Griq3VITwJZfA7bfD0Uf7ck45RAEsUltz5niXwzXX+LjV11+Hrl1jV5WbCgp8NMm99/p+7t8fnn46dlUpowAWqan1631Ws759fR7bhx6CMWOgZcvYleW+s86CyZP9H93/+3/ePzxrVuyq6k0BLLIjq1f7sKjdd/d5bL/3PZg+3cf7anazhtO3L0ycCLfcAm+84bdPP91/FllKASxSndJS+Pe/4Ywz/Cytq66C/ff304sffVRdDrE0aeJjhWfNgosv9j7ivfaCQw/15e5XroxdYa1YyMMp4GRrAwcODJMmTYpdRjxlZX4ixZtv+miGF1+EFSu8e+GHP/ThUAMGxK5StlRc7LOojR4Nn3zifcaHHgpDh/r1gAGw005pe3kzmxxCGFjn71cAC+RBAJeUeOuouBiWLPGB/vPmwezZHrxTp8KGDb5tp05wxBHe13jUUWn9A5YUCQHefReeegr++U9f/gm8i6h3bx+hssceUFTkZyp27uwLo7Zr5z/fOnYlKYAlJQZ26BAmDR+e/hfa1u9bxf1VHw9h60t5uV+XlVVeSkv9smkTfPWVXzZu9Mu6dbB2rX9dna5doU8f6NfPW0uDBvkfrPp2s9vy5T4vx5QpfsB05kz/Z1tSsvW2hYXQqpVfWrSA5s19CFzTpt7l0bixb1NQUHnZbTe49VYFsKTGwCZNwqSGWip9W+FWcX/Vx802vzRq5NdV/xgKCvyPpHHjyj+c5s390qqVdyO0aeOXDh2gY0cP3q5dfVvJD+Xl/unn88/9eulS/1S0cqX/k1671j8FbdxY+Y+8pMT/sZeWVv7DLy/31a1ffLHeAay5IMTtvTfkcheESKNG3vXQuXPsSr6mURAiIpEogEVEIlEfsABgZmuBmbHr2IH2wLLYRexANtQI2VFnNtTYJ4TQqq7frD5gqTCzPgcTGoKZTVKNqZENdWZLjfX5fnVBiIhEogAWEYlEASwV7o1dQA2oxtTJhjpzvkYdhBMRiUQtYBGRSBTAIiKRKIDznJkNM7OZZjbbzH4Ru54KZtbNzMab2XQzm2ZmFyf3tzOzl8zsk+S6bQbUWmBmU8xsXHJ7dzObmOzTx8ysSeT62pjZWDP72MxmmNmQTNuPZnZp8nOeamZjzKxZJuxHM7vfzJaa2dQq91W778z9Ian3QzPb4fylCuA8ZmYFwJ3Ad4G+wClmlinL+pYCl4cQ+gKDgfOT2n4BvBJC2AN4Jbkd28XAjCq3bwJuCyH0AlYCZ0apqtIdwAshhD2BffBaM2Y/mlkX4CJgYAihH1AAnExm7Me/AMO2uG9b++67wB7JZSRw9w6fPYSgS55egCHAi1VuXwVcFbuubdT6DPAd/Gy9zsl9nfETSGLW1TX5I/w2MA4w/Oytwur2cYT6WgNzSQ64V7k/Y/Yj0AWYD7TDTw4bBxyVKfsRKAKm7mjfAX8CTqluu21d1ALObxW/+BUWJPdlFDMrAvYDJgIdQwiLkocWAx0jlVXhduAKoDy5vQuwKoRQmtyOvU93B4qBB5JukvvMbCcyaD+GEBYCtwCfA4uA1cBkMms/VrWtfVfrvycFsGQ0M2sJPAlcEkJYU/Wx4M2MaOMozexYYGkIYXKsGmqgEBgA3B1C2A9YzxbdDRmwH9sCw/F/FrsBO7H1x/6MVN99pwDObwuBblVud03uywhm1hgP30dCCH9P7l5iZp2TxzsDS2PVBxwE/LeZfQb8De+GuANoY2YV86zE3qcLgAUhhInJ7bF4IGfSfjwCmBtCKA4hbAL+ju/bTNqPVW1r39X670kBnN/eBfZIjjY3wQ98PBu5JsCPKAOjgRkhhFurPPQsMCL5egTeNxxFCOGqEELXEEIRvu/+HUI4FRgPfD/ZLHaNi4H5ZtYnuWsoMJ0M2o9418NgM2uR/NwrasyY/biFbe27Z4HTk9EQg4HVVboqqher412XzLgARwOzgE+Ba2LXU6Wug/GPdh8C7yeXo/E+1leAT4CXgXaxa03qPQwYl3zdE3gHmA08ATSNXNu+wKRkXz4NtM20/QhcB3wMTAX+CjTNhP0IjMH7pTfhnybO3Na+ww/A3pn8LX2Ej+rY7vPX61RkMxuGf+QqAO4LIdxY5ycTEckzdQ7gZAzpLHxo0AL84+wpIYTpqStPRCR31WdC9gOA2SGEOQBm9jf8SOY2A7h9+/ahqKioHi8p9bVqlS/+2q3b5vdPnjx5WQihQ8Xt7zQ6UbM0iVTxUvkT21jOu+7qE8DVjXkbtOVGZjYSPyuE7t27M0kr70Z1xRUwatTWCyCb2bw4FYnkr7SPgggh3BtCGBhCGNihQ4cdf4Ok1aZN0Lhx7CpEBOoXwBk9hlSqV1ICTaJODSMiFeoTwBk7hlS2bcMG2Gmn2FWICNSjDziEUGpmFwAv4sPQ7g8hTEtZZZIWa9ZAy5axqxARqOey9CGE54DnUlSLNIDly2GXXWJXISKgU5HzzsKF0KlT7CpEBBTAeaWsDD77DL7xjdiVpFZBr90p6LV77DJEak0BnEemToXSUujXL3YlIgL17AOW7PLaa3596KFx60i1stlzY5cgUidqAeeRhx+Gvn23Pg1ZROJQAOeJCRPg3Xfh3HNjVxJHQZvWFLRpHbsMkc0ogPNASYkHb8eOcPrpsasRkQrqA85xIcA118D778PTT8POO8euKA5r1iz5avX2t2va1L8o98ngwqaSNFYl+U4t4BwWgs9+dsst3gIePjx2RSJSlVrAOWr1arjoInjoITj/fPjDH2JXFFfp4iU12i589RUAhbv3AKCsvX9ksGmf+nUrP4+7bEnMNSwlV6gFnIOeftpHOzz8MPzqVz7/byP9pEUyjlrAOeTtt+G3v4V//AP23tuD+Fvfil1Vdiqdm8xPXzHEOOkb/nI/bxm3mJpMqlxW5tsvWtyQ5UmOULsoy5WVwVNPwUEHwYEHwn/+Azfc4CteKHxFMptawFlqxgx47DF45BGYPRuKiryf98c/1nST6VDRN9z8Xe8LpnlzANYN7A5AeWNvGTdZXQpA089XAlD2yZyGLFOyjAI4i8ya5aH7+OM+r4OZn1Z8/fXwve9BoX6aIllFf7IZbN06eOMNePlleOkl+OgjD92DD/YDayecAJ07x64yv5QtX7HZ7WYLfBWuwh5+fveXvXYFYNX+fv3VtzsC0KK4HICd3/fRGKVzPkt7rZL5FMAZZNMmmDjRA/eVV/z04dJSX8PtoIPg9tvh+9+HLl1iVyoiqaAAjmjpUg/cisvbb8P69d7K3X9/uPxyOOIID9+ky1EyVOm8+QAUJtdtkxbx0qFdAVjZuwCAFXvuBkDjtX6983zvM97pJV/Nq3z9+gaqWDKBAriBfPWVnw48YYKH7YQJMDcZ4lRQAP37w4gRMHQoHHYYtGsXtVwRaQAK4DRYtw4+/NADd8oUv/7wQ58UB7wLYfBgPz140CBv7Wql4txS0SJud79fh4P2BaB4vxYArO/qc02s/qZvX3DQ3gC0nG8AdJyw1h9456MGqVfiUADX05IlmwftlCnwySc+DwN4S3a//eDiiz1sBw2Crl3j1iwimUEBXEMbN8L06T4SoeplcZUToIqKPGxPPRX23de/7trV+3Qlv9mb7wOw65t+u2SYnyXzxUH+J7ipm3882tTzSwBmH9AEgEZfDAGgw2T/j77zk5MACKWlDVC1pJsCeAtlZTBnztZBO3s2lPtIIpo187kWjjqqMmj32QfatIlbu4hkl7wO4KVLtw7aadNgwwZ/3MxXEO7fH04+2a/794devfzAmUhdNXnhXQCK/uW/SKtOOwCA4gP8T7JtDz+Trku3Rf74Pj4M5pPj+gPQ8l2/3fXvnwNQOn9BQ5QtKZYXAbxhgwfrlmG7tMqMgh06eLiedVZl0O61lw6OiUj65FwAr1jhB8Lee6/yetasyoNizZr5suzHHFMZtP37+3I9Ig2u3GdTa/PQ2wC0f83HD39+ol9/MshnYRvc7TMAju/ygT/+TR+n+NaRuwOwcpr3Ffccu86fV6MnskLWBnAI8MUXW4ft559XbtOtGwwYsHn3wTe+oe4DEckMWRPA5eXejfD66z4/wuuvw6JFlY/37g1DhvjqD/vt55f27ePVK1IXFeOHd7slGT88ZB8A/jO8HwDLh3if2FEd/My5YX28pftlb5+f+OEhgwGYPM37lIue8iPHTV6clPbapfYyNoA3bYLJkyvD9s03YaUfl6BLFz9bbPBgb+Husw+0ahW1XBGRWsuoAA7BJxR/6CGfcnHNGr+/d2+fbvGQQ3z6xaIija2V/GBve59vT+8iZsUJgwC45SjvIz5mvw8BGNnhNQB+1/1pADZ08362Pw08FIBxpw4AoMtT3lJu8dTEdJcuNbDDADazbsBDQEcgAPeGEO4ws3bAY0AR8BlwUghhZV2KmDMH/vpXD945c3zkwQknwHHH+dSLnTrV5VlFRDKbhYrhAdvawKwz0DmE8J6ZtQImA8cDZwArQgg3mtkvgLYhhCu391wDBw4MkyZt3hf1+9/Dz37mLdqhQ+H00721q+FfDcvMJocQBlbc/k6jE7f/iyEZYdXpPvph/fH+cfFHe7wDwMg23nJu0chbvCvKfEWPe1Z6C3rMdP9Rd3yimW/3d7WId+Sl8idS/rl7hy3gEMIiYFHy9VozmwF0AYYDhyWbPQi8Cmw3gLd0770eviecALfd5qMWRETyxQ5bwJttbFYEvA70Az4PIbRJ7jdgZcXtbanaAp4714eE7b+/H2Br0qRub0BSQy3g3FB8jreIC49bBsBPe/rkEz/e2UdVNDbvG15Z5qd73rjsQACenL4fAF3+5i3mZv94p4Eqzh7paAHXeFVkM2sJPAlcEkJYU/Wx4Cle7R+smY00s0lmNqm4uPjr+7t2hcMP9/G7zz1Xt+JFRLJZjVrAZtYYGAe8GEK4NblvJnBYCGFR0k/8agihz/aeZ8s+4LVr4TvfgXff9YltRoyA4cP9bDVpWGoB55hkmNCSC7xFvNMxPm3fWUX/AeBHrfz2V8FnVVtW7rOx3bjkCACef8/nnOj5mI8jLvz35IaoOqNFaQEn3QujgRkV4Zt4FhiRfD0CeKa2L96qFbzwAlx5pc/NcPLJPuJh5EjvlqiYfUxEJBfVZBTEwcAbwEdARSReDUwEHge6A/PwYWgrqn2SRHWjICqUlcGrr8KDD8KTT/oEOm3b+jC0Qw/1McADBkDjxrV5e1JTagHnuEbe97v0PB8F0fRon4nqpO7vAXBa6+SMuiQPPtnUGoBRC7xFPG1CTwB6jVkNQPn70xui6owSaxTEf4BtvfDQVBVSUODD0IYOhbvugmeegfHj/Uy4f/zDt2nRwk83rjgh44ADNFxNRLJXrUZB1Nf2WsDbs3ixnyFXMQ/EBx/4WXNm0KePz/swYEDlHBBa0LL21ALOT8vP9D7ijUf7cfVhRTMAOK6Nr+BRErzlPGVjEQCPf+ajJUpe84lWuj8yB4DSRVWWhslRUVrAmaBTJ/j+9/3qpOywAAALFklEQVQCsGoVvPWWH7ybMsXDecyYyu179KgM5AEDfNWK3XbT6csiklmyogVcE8uWeRhXXN57b/PFMdu2rZyScu+9/bpfP03iU0EtYAFYd6L3EX9xpB/u2WsPX2ljcLu5ADRKRptOWOnzEE+d2wWATi/4wZlWf5vQcMU2sHS0gHMmgKuzdq13V3zwgY+y+PBDmDrV769QVLT5xOz9+/vkP/l2sE8BLNUp/fb+AMz/jp8p1eybqwDo3savWxT68LUPFnoQh0/9oEzROD/Rw976oOGKTbO87YKoq1atfBTFwQdX3hcCzJu39fJEzz8PFQvNNmkCe+65dTBrhWMRSaWcbgHXxldfwcyZ3kquGswLqqx12KaNd1ts2ZWx887x6k4VtYClJhrtvScAXwz1I91r+nsLuGlLn+yncWNfYmn9fO/bazutEZ2eSQ7ULV7SoLWmmlrAadS0qQfq3ntvfv/Kld5tUTWUH3mkcq5i8FWSK5anr7ju1EmtZRHZPgXwDrRt6+OODzmk8r4QYP58by2//75f3nsPxo6t3GbXXSsDed99fdKhXr0UypLdyj/8GIBOPg88XTv5arYrDveDcsv2SX7BW/lBvJX9yvlyl28A0PGdrgA0flmnNVdQANeBGXTv7pdjj628f/VqD+UpUzyUp0yBW2/15ZXAxycPGuSXwYP9RJK2beO8BxGJT33AaVZSAtOnw6RJMHEiTJjgi4tW7PbevT2MBw2CAw/0LpBGNZ6jLnXUByypVDLsWwAs79uYr9r6r1LhBm8dN1vutzu+6TMXlE2bGaHC2lMfcBZq0qSyG+KnP/X71qzxQJ4wwUP5hRd8OSbwlZwPPxyOOMJPy+7ZU90WIrlKLeAMUDE07vXX4ZVX/LJwoT9WVFQ5R8awYenrslALWNKloG9vANb28V/e9Z38I16hDxWm1XwfSVH4qp/+THlZwxZYQ2oB5ygzD9qiIl8TLwQfElcRxk8+CaNH+8khRx4JJ53k8ya3bh27chGpD7WAs0BZmc978eST8Pjj8Pnn3rUxbJiH8Qkn1H8Se7WApaEU9vDFHzd805c7/7KtT/hTUOK/cq2n+Hjh0jmfNXxx2xF1SSKJp6DAD9T97nfw2Wfw9ttw/vkweTKcdppPPvSb38Dy5bErFZHaUAs4i5WX+5zJv/+9n0rdvDn85Cfw8597KNeGWsASS8EuflZdyd5FAJQXJKMl5vvk72UzZ0epa0tqActmGjXyg3PPPVe5pNO990LfvnD77d51ISKZSwGcI/r1g/vv9yk4DzsMLr3UxxVPmxa7MpHtK1u+grLlKygY/x4F49+j6cRZNJ046+vHC/r2pqBvbxq1akWjHJs/VgGcY3r0gHHj4NFHYc4cD+E334xdlYhUR8PQcpAZnHKKT8N5xBE+dO3ZZ727QiTTlVdM2D1z7Wb3V4yesI6+HFLZ7LkNWlc6qAWcw7p185M7evb05ZwqTu4QkcygAM5xHTvCU0/5fMdnnVU5B4VItimdN5/SefOxL0uwL0so6NCBgg4dYpdVLwrgPNCrF9xwgw9VGz8+djUiUkEBnCfOOcenw7zrrtiViNRP6YKFlC5YSFlxMWXFxbHLqRcFcJ5o1gxGjIBnnoH162NXIyKgAM4rRx7pC49OyN2Vw0WyigI4jwwZ4tfvvBO3DhFxCuA80rq1j4r49NPYlYikTzadMacAzjNFRT6dpYjEpzPh8kz79rBoUewqRNLn6zPpsoBawHmmbVtYuTJ2FSICtQhgMyswsylmNi65vbuZTTSz2Wb2mJk1SV+ZkiotWsDGjbGrEBGoXQv4YmBGlds3AbeFEHoBK4EzU1mYpEezZgpgkUxRowA2s67AMcB9yW0Dvg2MTTZ5EDg+HQVKajVuDJs2xa5CRKDmLeDbgSuA8uT2LsCqEEJpcnsB0KW6bzSzkWY2ycwmFWf5aYO5oLDQT8YQkfh2GMBmdiywNIQwuS4vEEK4N4QwMIQwsEOWz1yUCxo18rXkRCS+mgxDOwj4bzM7GmgG7AzcAbQxs8KkFdwV0GyzWUABLJI5dtgCDiFcFULoGkIoAk4G/h1COBUYD3w/2WwE8EzaqpSUMdOcwCKZoj7jgK8ELjOz2Xif8OjUlCTppAAWyRy1OhMuhPAq8Gry9RzggNSXJOlkFrsCEamgM+FERCJRAOcZtYBFMocCWEQkEgWwiEgkCmARkUgUwCIikSiARUQiUQCLiESiABYRiUQBLCISiQJYRCQSBbCISCQKYBGRSBTAIiKRKIBFRCJRAIuIRKIAFhGJRAEsIhKJAlhEJBIFsIhIJApgEZFIFMAiIpEogEVEIlEAi4hEogAWEYlEASwiEokCWEQkEgWwiEgkCmARkUgUwCIikSiARUQiqVEAm1kbMxtrZh+b2QwzG2Jm7czsJTP7JLlum+5iRURySU1bwHcAL4QQ9gT2AWYAvwBeCSHsAbyS3BYRkRraYQCbWWvgUGA0QAihJISwChgOPJhs9iBwfLqKFBHJRTVpAe8OFAMPmNkUM7vPzHYCOoYQFiXbLAY6VvfNZjbSzCaZ2aTi4uLUVC0ikgNqEsCFwADg7hDCfsB6tuhuCCEEIFT3zSGEe0MIA0MIAzt06FDfekVEckZNAngBsCCEMDG5PRYP5CVm1hkguV6anhJFRHLTDgM4hLAYmG9mfZK7hgLTgWeBEcl9I4Bn0lKhiEiOKqzhdhcCj5hZE2AO8GM8vB83szOBecBJ6SlRRCQ31SiAQwjvAwOreWhoassREckfOhNORCQSBbCISCQKYBGRSBTAIiKRKIBFRCJRAIuIRKIAFhGJRAEsIhKJAlhEJBIFsIhIJApgEZFIFMAiIpEogEVEIlEAi4hEogAWEYlEASwiEokCWEQkEgWwiEgkCmARkUgUwCIikSiARUQiUQCLiESiABYRiUQBLCISiQJYRCQSBbCISCQKYBGRSBTAIiKRKIBFRCJRAIuIRFKjADazS81smplNNbMxZtbMzHY3s4lmNtvMHjOzJukuVkQkl+wwgM2sC3ARMDCE0A8oAE4GbgJuCyH0AlYCZ6azUBGRXFPTLohCoLmZFQItgEXAt4GxyeMPAsenvjwRkdy1wwAOISwEbgE+x4N3NTAZWBVCKE02WwB0qe77zWykmU0ys0nFxcWpqVpEJAfUpAuiLTAc2B3YDdgJGFbTFwgh3BtCGBhCGNihQ4c6Fyoikmtq0gVxBDA3hFAcQtgE/B04CGiTdEkAdAUWpqlGEZGcVJMA/hwYbGYtzMyAocB0YDzw/WSbEcAz6SlRRCQ31aQPeCJ+sO094KPke+4FrgQuM7PZwC7A6DTWKSKScwp3vAmEEH4N/HqLu+cAB6S8IhGRPKEz4UREIlEAi4hEogAWEYlEASwiEokCWEQkEgWwiEgkCmARkUgUwCIikSiARUQiUQCLiESiABYRiUQBLCISiQJYRCQSBbCISCQKYBGRSBTAIiKRKIBFRCJRAIuIRKIAFhGJRAEsIhKJAlhEJBIFsIhIJApgEZFIFMAiIpEogEVEIlEAi4hEogAWEYlEASwiEokCWEQkEgWwiEgkCmARkUgUwCIikSiA88zgwXDxxbGrEBEACyE03IuZFQPzGuwFpTZ6hBA6xC5CJJ80aACLiEgldUGIiESiABYRiUQBLCISiQJYRCQSBbCISCQKYBGRSBTAIiKRKIBFRCJRAIuIRKIAFhGJRAEsIhKJAlhEJBIFsIhIJApgEZFIFMAiIpEogEVEIlEAi4hEogAWEYlEASwiEokCWEQkEgWwiEgkCmARkUj+P3OwqFnETUiQAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 360x360 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#%% Smooth OT with KL regularization\n",
    "\n",
    "lambd = 2e-3\n",
    "Gsm = ot.smooth.smooth_ot_dual(a, b, M, lambd, reg_type='kl')\n",
    "\n",
    "pl.figure(5, figsize=(5, 5))\n",
    "ot.plot.plot1D_mat(a, b, Gsm, 'OT matrix Smooth OT KL reg.')\n",
    "\n",
    "pl.show()\n",
    "\n",
    "\n",
    "#%% Smooth OT with KL regularization\n",
    "\n",
    "lambd = 1e-1\n",
    "Gsm = ot.smooth.smooth_ot_dual(a, b, M, lambd, reg_type='l2')\n",
    "\n",
    "pl.figure(6, figsize=(5, 5))\n",
    "ot.plot.plot1D_mat(a, b, Gsm, 'OT matrix Smooth OT l2 reg.')\n",
    "\n",
    "pl.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
